AliPhysics  master (3d17d9d)
AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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  * appeuear 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 /* $Id$ */
17 
18 //
19 //
20 // XicPlus2XiPiPi analysis code
21 //
22 // Input: AOD
23 // Output: TTree or THnSparse (mass vs pT vs Centrality)
24 //
25 // Cuts:
26 // TTree: very loose cut
27 // THnSparse: One THnSparse is created per cut. One cut is specified by
28 // an array of bits, each bit corresponds to a cut in "Cut" function.
29 //
30 //-------------------------------------------------------------------------
31 //
32 // Authors: Y.S Watanabe(a), G. Luparello (b)
33 //
34 // (a) CNS, the University of Tokyo
35 // (b) INFN, Trieste
36 //
37 // Contatcs: wyosuke@cns.s.u-tokyo.ac.jp, gluparel@cern.ch
38 //-------------------------------------------------------------------------
39 
40 #include <TSystem.h>
41 #include <TParticle.h>
42 #include <TParticlePDG.h>
43 #include <THnSparse.h>
44 #include <TH1F.h>
45 #include <TH2F.h>
46 #include <TH3F.h>
47 #include <TLorentzVector.h>
48 #include <TTree.h>
49 #include "TROOT.h"
50 #include <TDatabasePDG.h>
51 #include <AliAnalysisDataSlot.h>
52 #include <AliAnalysisDataContainer.h>
53 #include "AliMCEvent.h"
54 #include "AliAnalysisManager.h"
55 #include "AliAODMCHeader.h"
56 #include "AliAODHandler.h"
57 #include "AliLog.h"
58 #include "AliExternalTrackParam.h"
59 #include "AliAODVertex.h"
60 #include "AliESDVertex.h"
61 #include "AliAODRecoDecay.h"
62 #include "AliAODRecoDecayHF.h"
64 #include "AliESDtrack.h"
65 #include "AliAODTrack.h"
66 #include "AliAODv0.h"
67 #include "AliAODcascade.h"
68 #include "AliAODMCParticle.h"
69 #include "AliAnalysisTaskSE.h"
71 #include "AliPIDResponse.h"
72 #include "AliPIDCombined.h"
73 #include "AliAODPidHF.h"
74 #include "AliInputEventHandler.h"
75 #include "AliESDtrackCuts.h"
76 #include "AliNeutralTrackParam.h"
77 #include "AliKFParticle.h"
78 #include "AliKFVertex.h"
79 #include "AliExternalTrackParam.h"
80 #include "AliCentrality.h"
81 #include "AliVertexerTracks.h"
83 #include "AliVertexingHFUtils.h"
84 
85 using std::cout;
86 using std::endl;
87 
91 
92 //__________________________________________________________________________
95  fUseMCInfo(kFALSE),
96  fFillSignalOnly(kFALSE),
97  fFillBkgOnly(kFALSE),
98  fOutput(0),
99  fOutputAll(0),
100  fListCuts(0),
101  fCEvents(0),
102  fHTrigger(0),
103  fHCentrality(0),
104  fAnalCuts(0),
105  fIsEventSelected(kFALSE),
106  fWriteVariableTree(kFALSE),
107  fFillSparse(kFALSE),
108  fVariablesTree(0),
109  fReconstructPrimVert(kFALSE),
110  fIsMB(kFALSE),
111  fIsSemi(kFALSE),
112  fIsCent(kFALSE),
113  fIsINT7(kFALSE),
114  fIsEMC7(kFALSE),
115  fCandidateVariables(),
116  fVtx1(0),
117  fV1(0),
118  fBzkG(0),
119  fCentrality(0),
120  //fTriggerCheck(0),
121  fHistoXicMass(0x0),
122  fSparseXicMass(0x0),
123  fHistoMCSpectrumAccXic(0),
124  fHistoDcaPi1Pi2(0),
125  fHistoDcaPi1Casc(0),
126  fHistoDcaPi2Casc(0),
127  fHistoLikeDecayLength(0),
128  fHistoLikeDecayLengthXY(0),
129  fHistoXicCosPAXY(0),
130  fHistoXiMass(0),
131  fHistoCascDcaXiDaughters(0),
132  fHistoCascDcaV0Daughters(0),
133  fHistoCascDcaV0ToPrimVertex(0),
134  fHistoCascDcaPosToPrimVertex(0),
135  fHistoCascDcaNegToPrimVertex(0),
136  fHistoCascDcaBachToPrimVertex(0),
137  fHistoCascCosPAXiPrim(0),
138  fHistoXiPt(0),
139  fHistoPiPt(0),
140  fHistoPid0(0),
141  fHistonSigmaTPCpi(0),
142  fHistonSigmaTOFpi(0),
143  fHistoProbPion(0),
144  fHistoXiMassvsPtRef1(0),
145  fHistoXiMassvsPtRef2(0),
146  fHistoXiMassvsPtRef3(0),
147  fHistoXiMassvsPtRef4(0),
148  fHistoXiMassvsPtRef5(0),
149  fHistoXiMassvsPtRef6(0),
150  fHistoPiPtRef(0),
151  fHistoPiEtaRef(0),
152  fQAHistoNSelectedTracks(0),
153  fQAHistoNSelectedCasc(0),
154  fQAHistoDCApi1pi2(0),
155  fQAHistoAODPrimVertX(0),
156  fQAHistoAODPrimVertY(0),
157  fQAHistoAODPrimVertZ(0),
158  fQAHistoRecoPrimVertX(0),
159  fQAHistoRecoPrimVertY(0),
160  fQAHistoRecoPrimVertZ(0),
161  fQAHistoSecondaryVertexX(0),
162  fQAHistoSecondaryVertexY(0),
163  fQAHistoSecondaryVertexZ(0),
164  fQAHistoSecondaryVertexXY(0),
165  fCounter(0)
166 {
167  //
168  // Default Constructor.
169  //
170 }
171 
172 //___________________________________________________________________________
175  Bool_t writeVariableTree, Bool_t fillSparse) :
176  AliAnalysisTaskSE(name),
177  fUseMCInfo(kFALSE),
178  fFillSignalOnly(kFALSE),
179  fFillBkgOnly(kFALSE),
180  fOutput(0),
181  fOutputAll(0),
182  fListCuts(0),
183  fCEvents(0),
184  fHTrigger(0),
185  fHCentrality(0),
186  fAnalCuts(analCuts),
187  fIsEventSelected(kFALSE),
188  fWriteVariableTree(writeVariableTree),
189  fFillSparse(fillSparse),
190  fVariablesTree(0),
191  fReconstructPrimVert(kFALSE),
192  fIsMB(kFALSE),
193  fIsSemi(kFALSE),
194  fIsCent(kFALSE),
195  fIsINT7(kFALSE),
196  fIsEMC7(kFALSE),
198  fVtx1(0),
199  fV1(0),
200  fBzkG(0),
201  fCentrality(0),
202  //fTriggerCheck(0),
203  fHistoXicMass(0x0),
204  fSparseXicMass(0x0),
206  fHistoDcaPi1Pi2(0),
207  fHistoDcaPi1Casc(0),
208  fHistoDcaPi2Casc(0),
211  fHistoXicCosPAXY(0),
212  fHistoXiMass(0),
220  fHistoXiPt(0),
221  fHistoPiPt(0),
222  fHistoPid0(0),
225  fHistoProbPion(0),
232  fHistoPiPtRef(0),
233  fHistoPiEtaRef(0),
247  fCounter(0)
248 {
249  //
250  // Constructor. Initialization of Inputs and Outputs
251  //
252  Info("AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks","Calling Constructor");
253 
254  DefineOutput(1,TList::Class());
255  DefineOutput(2,TList::Class());
256  if(writeVariableTree){
257  DefineOutput(3,TTree::Class());
258  }else{
259  DefineOutput(3,TList::Class());
260  }
261  DefineOutput(4,AliNormalizationCounter::Class());
262 }
263 
264 //___________________________________________________________________________
266  //
267  // destructor
268  //
269  Info("~AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks","Calling Destructor");
270 
271 
272  if (fOutput) {
273  delete fOutput;
274  fOutput = 0;
275  }
276 
277  if (fOutputAll) {
278  delete fOutputAll;
279  fOutputAll = 0;
280  }
281 
282  if (fListCuts) {
283  delete fListCuts;
284  fListCuts = 0;
285  }
286 
287  if (fAnalCuts) {
288  delete fAnalCuts;
289  fAnalCuts = 0;
290  }
291 
292  if (fVariablesTree) {
293  delete fVariablesTree;
294  fVariablesTree = 0;
295  }
296 
298  delete fCandidateVariables;
300  }
301 
302  if(fHistoXicMass) delete fHistoXicMass;
303  if(fSparseXicMass) delete fSparseXicMass;
304 
306 
307  if(fCounter){
308  delete fCounter;
309  fCounter=0;
310  }
311 }
312 
313 //_________________________________________________
315  //
316  // Initialization
317  //
318  //
319 
320  //Copied from $ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF.C
321 
322  fIsEventSelected=kFALSE;
323 
324  if (fDebug > 1) AliInfo("Init");
325 
326  fListCuts = new TList();
327  fListCuts->SetOwner();
328  fListCuts->SetName("ListCuts");
330  PostData(2,fListCuts);
331 
332  return;
333 }
334 
335 //_________________________________________________
337 {
338  //
339  // UserExec code
340  //
341 
342  if (!fInputEvent) {
343  AliError("NO EVENT FOUND!");
344  return;
345  }
346  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
347 
348  fCEvents->Fill(1);
349  //------------------------------------------------
350  // First check if the event has proper vertex and B
351  //------------------------------------------------
352  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
353  if (!fVtx1) return;
354 
355  Double_t pos[3],cov[6];
356  fVtx1->GetXYZ(pos);
357  fVtx1->GetCovarianceMatrix(cov);
358  fQAHistoAODPrimVertX->Fill(pos[0]);
359  fQAHistoAODPrimVertY->Fill(pos[1]);
360  fQAHistoAODPrimVertZ->Fill(pos[2]);
361 
362  fV1 = new AliESDVertex(pos,cov,100.,100,fVtx1->GetName());
363 
364  fBzkG = (Double_t)aodEvent->GetMagneticField();
365  AliKFParticle::SetField(fBzkG);
366  if (TMath::Abs(fBzkG)<0.001) {
367  delete fV1;
368  return;
369  }
370  fCEvents->Fill(2);
371 
373 
374  //------------------------------------------------
375  // Event selection
376  //------------------------------------------------
377  Bool_t fIsTriggerNotOK = fAnalCuts->IsEventRejectedDueToTrigger();
378  if(!fIsTriggerNotOK) fCEvents->Fill(3);
379 
381  if(!fIsEventSelected) {
382  //cout<<"Why: "<<fAnalCuts->GetWhyRejection()<<endl;
383  delete fV1;
384  return;
385  }
386 
387  //cout<<fabs(aodEvent->GetPrimaryVertex()->GetZ()-aodEvent->GetPrimaryVertexSPD()->GetZ())<<endl;
388 
389  fCEvents->Fill(4);
390 
391  fIsMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB)==(AliVEvent::kMB);
392  fIsSemi=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kSemiCentral)==(AliVEvent::kSemiCentral);
393  fIsCent=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kCentral)==(AliVEvent::kCentral);
394  fIsINT7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kINT7)==(AliVEvent::kINT7);
395  fIsEMC7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kEMC7)==(AliVEvent::kEMC7);
396  //this trigger check is not used: commentig out for the moment everywhere in the task
397  // fTriggerCheck = fIsMB+2*fIsSemi+4*fIsCent+8*fIsINT7+16*fIsEMC7;
398  if(fIsMB) fHTrigger->Fill(1);
399  if(fIsSemi) fHTrigger->Fill(2);
400  if(fIsCent) fHTrigger->Fill(3);
401  if(fIsINT7) fHTrigger->Fill(4);
402  if(fIsEMC7) fHTrigger->Fill(5);
403  if(fIsMB|fIsSemi|fIsCent) fHTrigger->Fill(7);
404  if(fIsINT7|fIsEMC7) fHTrigger->Fill(8);
405  if(fIsMB&fIsSemi) fHTrigger->Fill(10);
406  if(fIsMB&fIsCent) fHTrigger->Fill(11);
407  if(fIsINT7&fIsEMC7) fHTrigger->Fill(12);
408 
409  AliCentrality *cent = aodEvent->GetCentrality();
410  fCentrality = cent->GetCentralityPercentile("V0M");
411  fHCentrality->Fill(fCentrality);
412 
413  //------------------------------------------------
414  // MC analysis setting
415  //------------------------------------------------
416  TClonesArray *mcArray = 0;
417  AliAODMCHeader *mcHeader=0;
418  if (fUseMCInfo) {
419  // MC array need for maching
420  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
421  if (!mcArray) {
422  AliError("Could not find Monte-Carlo in AOD");
423  delete fV1;
424  return;
425  }
426  fCEvents->Fill(6); // in case of MC events
427 
428  // load MC header
429  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
430  if (!mcHeader) {
431  AliError("AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::UserExec: MC header branch not found!\n");
432  delete fV1;
433  return;
434  }
435  fCEvents->Fill(7); // in case of MC events
436 
437  Double_t zMCVertex = mcHeader->GetVtxZ();
438  if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()) {
439  AliDebug(2,Form("Event rejected: abs(zVtxMC)=%f > fAnalCuts->GetMaxVtxZ()=%f",zMCVertex,fAnalCuts->GetMaxVtxZ()));
440  delete fV1;
441  return;
442  } else {
443 
444  //DOUBT: isEventSelected before or after LoopOverGenParticles?
445 
446  LoopOverGenParticles(mcArray);
447  cout<<"Loop over gen particles done!!!"<<endl;
448  fCEvents->Fill(17); // in case of MC events
449  }
450  }
451 
452  //------------------------------------------------
453  // Check if the event has cascade candidate
454  //------------------------------------------------
455  Int_t ncasc = aodEvent->GetNumberOfCascades();
456  Int_t nselecasc = 0.;
457  for(Int_t ic=0;ic<ncasc;ic++){
458  AliAODcascade *casc = aodEvent->GetCascade(ic);
459  if(!fAnalCuts) continue;
460  if(fAnalCuts->SingleCascadeCuts(casc,pos )) nselecasc++;
461  }
462 
463  if(nselecasc==0){
464  delete fV1;
465  return;
466  }
467 
468  fCEvents->Fill(5); //counter of the events with selected cascades
469 
470  //------------------------------------------------
471  // Main analysis done in this function
472  //------------------------------------------------
473  MakeAnalysis(aodEvent, mcArray);
474 
475  PostData(1,fOutput);
476  if(fWriteVariableTree){
477  PostData(3,fVariablesTree);
478  }else{
479  PostData(3,fOutputAll);
480  }
481  PostData(4,fCounter);
482  fIsEventSelected=kFALSE;
483 
484  delete fV1;
485  return;
486 }
487 
488 //________________________________________ terminate ___________________________
490 {
491  // The Terminate() function is the last function to be called during
492  // a query. It always runs on the client, it can be used to present
493  // the results graphically or save the results to file.
494 
495  //AliInfo("Terminate","");
496  AliAnalysisTaskSE::Terminate();
497 
498  fOutput = dynamic_cast<TList*> (GetOutputData(1));
499  if (!fOutput) {
500  AliError("fOutput not available");
501  return;
502  }
503 
504  if(!fWriteVariableTree){
505  fOutputAll = dynamic_cast<TList*> (GetOutputData(3));
506  if (!fOutputAll) {
507  AliError("fOutputAll not available");
508  return;
509  }
510  }
511 
512  return;
513 }
514 
515 //___________________________________________________________________________
517 {
518  //
519  // UserCreateOutputObjects
520  //
521 
522  AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
523 
524  //------------------------------------------------
525  // output object setting
526  //------------------------------------------------
527  fOutput = new TList();
528  fOutput->SetOwner();
529  fOutput->SetName("chist0");
530  DefineGeneralHistograms(); // define general histograms
531  PostData(1,fOutput);
532 
533 
534  if (fWriteVariableTree) {
536  PostData(3,fVariablesTree);
537  }else{
538  fOutputAll = new TList();
539  fOutputAll->SetOwner();
540  fOutputAll->SetName("anahisto");
541  DefineAnalysisHistograms(); // define general histograms
542  PostData(3,fOutputAll);
543  }
544 
545  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(4)->GetContainer()->GetName()));
546  fCounter->Init();
547  PostData(4,fCounter);
548  return;
549 }
550 
551 //________________________________________________________________________
553 (
554  AliAODEvent *aodEvent, TClonesArray *mcArray
555  )
556 {
557  //
558  // Main analysis part called from UserExec
559  //
560 
561  Int_t nCascades= aodEvent->GetNumberOfCascades();
562  if (nCascades==0) {
563  return;
564  }
565  Int_t nTracks= aodEvent->GetNumberOfTracks();
566  if (nTracks==0) {
567  return;
568  }
569 
570  //------------------------------------------------
571  // Select good track before hand to save time
572  //------------------------------------------------
573  Bool_t seleTrkFlags[nTracks];
574  Int_t nSeleTrks=0;
575  SelectTrack(aodEvent,nTracks,nSeleTrks,seleTrkFlags); //select candidates pions
576  fQAHistoNSelectedTracks->Fill(nSeleTrks);
577 
578  Bool_t seleCascFlags[nCascades];
579  Int_t nSeleCasc=0;
580  SelectCascade(aodEvent,nCascades,nSeleCasc,seleCascFlags);
581  fQAHistoNSelectedCasc->Fill(nSeleCasc);
582 
583  Int_t usedmclab[20];//Used Xic Label: Assuming there are less than 20 Xic/evt
584  Int_t nusedmclab[20];//Number of times the Xic label is used: Assuming there are less than 20 Xic/evt
585  for(Int_t i=0;i<20;i++) {
586  usedmclab[i]=-9999;
587  nusedmclab[i]=0;
588  }
589 
590  //------------------------------------------------
591  // Cascade loop
592  //------------------------------------------------
593  for (Int_t icasc = 0; icasc<nCascades; icasc++) {
594  if(!seleCascFlags[icasc]) continue;
595  fCEvents->Fill(18);
596 
597  AliAODcascade *casc = aodEvent->GetCascade(icasc);
598  if(!casc) continue;
599 
600 
601  AliAODTrack *cptrack = (AliAODTrack*)(casc->GetDaughter(0));
602  AliAODTrack *cntrack = (AliAODTrack*)(casc->GetDaughter(1));
603  AliAODTrack *cbtrack = (AliAODTrack*)(casc->GetDecayVertexXi()->GetDaughter(0));
604  if(!cptrack || !cntrack || !cbtrack) continue;
605 
606  Int_t cpid = cptrack->GetID();
607  Int_t cnid = cntrack->GetID();
608  Int_t cbid = cbtrack->GetID();
609 
610  if(cptrack->Charge()==0) continue;
611  if(cntrack->Charge()==0) continue;
612  if(cbtrack->Charge()==0) continue;
613 
614  Short_t charge_casc = cptrack->Charge() + cntrack->Charge() + cbtrack->Charge();
615 
616  //------------------------------------------------
617  // Track1 loop
618  //------------------------------------------------
619  for (Int_t itrk1 = 0; itrk1<nTracks-1; itrk1++) {
620  if(!seleTrkFlags[itrk1]) continue;
621  AliAODTrack *trk1 = (AliAODTrack*)aodEvent->GetTrack(itrk1);
622  if(!trk1||trk1->GetID()<0) continue;
623 
624  //------------------------------------------------
625  // Track2 loop
626  //------------------------------------------------
627  for (Int_t itrk2 = itrk1+1; itrk2<nTracks; itrk2++) {
628  if(!seleTrkFlags[itrk2]) continue;
629  AliAODTrack *trk2 = (AliAODTrack*)aodEvent->GetTrack(itrk2);
630  if(!trk2||trk2->GetID()<0) continue;
631 
632  if(!SelectLikeSign(trk1,trk2)) continue;
633  fCEvents->Fill(19);
634 
635  Int_t lpid1 = trk1->GetID();
636  Int_t lpid2 = trk2->GetID();
637  if((cpid==lpid1)||(cpid==lpid2)||(cnid==lpid1)||(cnid==lpid2)||(cbid==lpid1)||(cbid==lpid2)) continue;
638 
639 
640  Short_t charge_like1 = trk1->Charge();
641  Short_t charge_like2 = trk2->Charge();
642  Bool_t ok_charge = kFALSE;
643  if((charge_casc==-1)&&(charge_like1==1)&&(charge_like2==1)) ok_charge = kTRUE;
644  if((charge_casc==1)&&(charge_like1==-1)&&(charge_like2==-1)) ok_charge = kTRUE;
645  if(!ok_charge) continue;
646  fCEvents->Fill(20);
647  //------------------------------------------------
648  // Roughly select good candidates to speed up
649  //------------------------------------------------
650  if(!fAnalCuts->SelectWithRoughCuts(casc,trk1,trk2)) continue;
651  fCEvents->Fill(21);
652  //------------------------------------------------
653  // Secondary vertex calculation
654  //------------------------------------------------
655  Double_t dispersion;
656  AliAODVertex *secVert = CallReconstructSecondaryVertex(trk1,trk2,dispersion);
657  if(!secVert) continue;
658  fCEvents->Fill(22);
659  fQAHistoSecondaryVertexX->Fill(secVert->GetX());
660  fQAHistoSecondaryVertexY->Fill(secVert->GetY());
661  fQAHistoSecondaryVertexZ->Fill(secVert->GetZ());
662  fQAHistoSecondaryVertexXY->Fill(secVert->GetX()*secVert->GetX()+secVert->GetY()*secVert->GetY());
663 
664  AliAODRecoCascadeHF3Prong *xicobj = MakeCascadeHF3Prong(casc,trk1,trk2,aodEvent,secVert,dispersion);
665  if(!xicobj) {
666  delete secVert;
667  continue;
668  }
669  fCEvents->Fill(23);
670 
671  AliAODMCParticle *mcxic = 0;
672  AliAODMCParticle *mcdaughter1 = 0;
673  AliAODMCParticle *mcdaughter2 = 0;
674  AliAODMCParticle *mcdaughterxi = 0;
675  AliAODMCParticle *mcdaughterLambda = 0;
676  AliAODMCParticle *mcdaughterPionFromLambda =0;
677  AliAODMCParticle *mcdaughterProtonFromLambda=0;
678  AliAODMCParticle *mcdaughterPionFromXi=0;
679  Int_t nFound=0;
680 
681  Int_t mclabxic = 0;
682  Int_t nmclabxic = 0;
683  Bool_t isXic = kFALSE;
684 
685  if(fUseMCInfo)
686  {
687  //Note: this can be enclosed in a MatchRecoCandToMC() function
688  //Input will be the AliAODRecoCascadeHF3Prong and mcArray and return the AODMCParticle
689 
690  Int_t pdgDg[3]={211,3312,211};
691  Int_t pdgDgcasc[2]={211,3122};
692  Int_t pdgDgv0[2]={2212,211};
693  mclabxic = xicobj->MatchToMC(4232,pdgDg[1],pdgDg,pdgDgcasc,pdgDgv0,mcArray);
694  if(mclabxic>-1){
695  mcxic = (AliAODMCParticle*) mcArray->At(mclabxic);
696  if (mcxic){
697  Int_t checkOrigin=AliVertexingHFUtils::CheckOrigin(mcArray,mcxic,kTRUE);
698  Bool_t isInAcc=kTRUE;
699  if(fAnalCuts){
700  if(!fAnalCuts->IsInFiducialAcceptance(mcxic->Pt(),mcxic->Y())) isInAcc=kFALSE;
701  }else{
702  if(TMath::Abs(mcxic->Y())>0.8) isInAcc=kFALSE;
703  }
704 
705  for(Int_t ia=0;ia<20;ia++){ //this is to check how many Xic there are per event.. I guess always only one. I am not sure that this is really needed
706  if(usedmclab[ia]==mclabxic){
707  nusedmclab[ia]++;
708  nmclabxic = nusedmclab[ia];
709  break;
710  }
711  if(usedmclab[ia]==-9999){
712  usedmclab[ia]=mclabxic;
713  nusedmclab[ia]++;
714  nmclabxic = nusedmclab[ia];
715  break;
716  }
717  }
718  Int_t pi_counter = 0;
719  for(Int_t idau=mcxic->GetDaughterFirst();idau<mcxic->GetDaughterLast()+1;idau++)
720  {
721  //cout<<idau<<endl;
722  if(idau<0) break;
723  AliAODMCParticle *mcpart = (AliAODMCParticle*) mcArray->At(idau);
724  Int_t pdgcode = TMath::Abs(mcpart->GetPdgCode());
725  if(pdgcode==211 && pi_counter==0){
726  mcdaughter1 = mcpart;
727  pi_counter=1;
728  nFound++;
729  }else if(pdgcode==211 && pi_counter==1){
730  mcdaughter2 = mcpart;
731  pi_counter=2;
732  nFound++;
733  }else if(pdgcode==3312){
734  mcdaughterxi = mcpart;
735  Int_t nXiDau= mcpart->GetNDaughters();
736  if(nXiDau!=2) break; //no correct Xi decay
737  Int_t indFirstXiDau=mcpart->GetDaughterLabel(0);
738  for(Int_t XiDau=0; XiDau<2; XiDau++){
739  Int_t indXiDau=indFirstXiDau+XiDau;
740  if(indXiDau<0) break;
741  AliAODMCParticle* Xidau=dynamic_cast<AliAODMCParticle*>(mcArray->At(indXiDau));
742  if(!Xidau) break;
743  Int_t pdgXidau=Xidau->GetPdgCode();
744  if(TMath::Abs(pdgXidau)==3122){
745  cout<< "-- I am a lambda" <<endl;
746  Int_t nLambdaDau=Xidau->GetNDaughters();
747  if (nLambdaDau!=2) break;
748  Int_t indFirstLambdaDau=Xidau->GetDaughterLabel(0);
749  for(Int_t LambdaDau=0; LambdaDau<2; LambdaDau++){
750  Int_t indLambdaDau=indFirstLambdaDau+LambdaDau;
751  AliAODMCParticle* Lambdadau=dynamic_cast<AliAODMCParticle*>(mcArray->At(indLambdaDau));
752  if(!Lambdadau) break;
753  Int_t pdgLambdadau=Lambdadau->GetPdgCode();
754  if(TMath::Abs(pdgLambdadau)==2212){
755  cout<< "-- I am a proton from Lambda" <<endl;
756  mcdaughterProtonFromLambda=Lambdadau;
757  nFound++;
758  } else if(TMath::Abs(pdgLambdadau)==211){
759  cout<< "-- I am a pion from Lambda" <<endl;
760  mcdaughterPionFromLambda=Lambdadau;
761  nFound++;
762  }
763  mcdaughterLambda=Xidau;
764  }
765  } else if(TMath::Abs(pdgXidau)==211) {
766  cout<< "-- I am a pion from Xi" <<endl;
767  mcdaughterPionFromXi=Xidau;
768  nFound++;
769  }
770  }
771  }
772  }
773  if(nFound==5){
774  if(isInAcc){
775  if ((TMath::Abs(mcdaughter1->Eta())>0.9) || (TMath::Abs(mcdaughter2->Eta())>0.9) || (TMath::Abs(mcdaughterProtonFromLambda->Eta())>0.9) || (TMath::Abs(mcdaughterPionFromLambda->Eta())>0.9) || (TMath::Abs(mcdaughterPionFromXi->Eta())>0.9)){ //to check also the y acceptance of Xi and Lambda?
776  isInAcc=kFALSE;
777  }
778  }
779  if(isInAcc){
780  fHistoMCSpectrumAccXic->Fill(mcxic->Pt(),kReco,checkOrigin);
781  isXic=kTRUE;
783  fHistoMCSpectrumAccXic->Fill(mcxic->Pt(),kRecoPID,checkOrigin);
784  }
785  fAnalCuts->SetUsePID(kFALSE);
787  fHistoMCSpectrumAccXic->Fill(mcxic->Pt(),kRecoCuts,checkOrigin);
788  }
789  fAnalCuts->SetUsePID(kTRUE);
790  }
791  }
792  }
793  //cout<<"pointer: "<<mcdaughter1<<" "<<mcdaughter2<<" "<<mcdaughterxi<<endl;
794  //cout<<"1: "<<mcdaughter1->GetPdgCode()<<endl;
795  //cout<<"2: "<<mcdaughter2->GetPdgCode()<<endl;
796  //cout<<"3: "<<mcdaughterxi->GetPdgCode()<<endl;
797  }
798  } //if useMCinfo
799 
800  if(!fAnalCuts->IsInFiducialAcceptance(xicobj->Pt(),xicobj->Y())) continue;
801 
802  FillROOTObjects(xicobj,mcxic,mcdaughter1,mcdaughter2,mcdaughterxi,nmclabxic,isXic); //AliAODRecoCascadeHF3Prong,
803  xicobj->GetSecondaryVtx()->RemoveDaughters();
804  xicobj->UnsetOwnPrimaryVtx();
805  delete xicobj;xicobj=NULL;
806  delete secVert;
807 
808  }//trk2
809  }//trk1
810  }//casc
811 }
812 
813 //________________________________________________________________________
814 void AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::FillROOTObjects(AliAODRecoCascadeHF3Prong *xicobj, AliAODMCParticle *mcpart, AliAODMCParticle *mcdaughter1, AliAODMCParticle *mcdaughter2, AliAODMCParticle *mcdaughterxi, Int_t mcnused, Bool_t isXic)
815 {
816  //
817  // Fill histogram or Tree depending on fWriteVariableTree flag
818  //
819 
820  AliAODTrack *part1 = xicobj->GetBachelor1();
821  AliAODTrack *part2 = xicobj->GetBachelor2();
822  AliAODcascade *casc = xicobj->GetCascade();
823 
824  Double_t dca[3];
825  xicobj->GetDCAs(dca);
826 
827  Double_t nSigmaTPCpi1=-9999.;
828  Double_t nSigmaTPCpi2=-9999.;
829  Double_t nSigmaTOFpi1=-9999.;
830  Double_t nSigmaTOFpi2=-9999.;
831  Double_t probPion1=-9999.;
832  Double_t probPion2=-9999.;
833 
834  if(fAnalCuts->GetIsUsePID()) {
835  nSigmaTPCpi1 = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTPC(part1,AliPID::kPion);
836  nSigmaTPCpi2 = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTPC(part2,AliPID::kPion);
837  nSigmaTOFpi1 = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTOF(part1,AliPID::kPion);
838  nSigmaTOFpi2 = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTOF(part2,AliPID::kPion);
839 
841  probPion1 = fAnalCuts->GetPionProbabilityTPCTOF(part1);
842  probPion2 = fAnalCuts->GetPionProbabilityTPCTOF(part2);
843  }
844  }
845 
846 
847  if(fWriteVariableTree){
848  if ((fUseMCInfo && fFillSignalOnly && !fFillBkgOnly && !isXic) || (fUseMCInfo && !fFillSignalOnly && fFillBkgOnly && isXic)) return;
849  else {
850 
851  // if (!fUseMCInfo || (fUseMCInfo && fFillSignalOnly && isXiC) || (fUseMCInfo && !fFillSignalOnly)){
852  fCandidateVariables[ 0] = xicobj->InvMassPiXiPi();
853  fCandidateVariables[ 1] = xicobj->Pt();
854  fCandidateVariables[ 2] = part1->Px();
855  fCandidateVariables[ 3] = part1->Py();
856  fCandidateVariables[ 4] = part1->Pz();
857  fCandidateVariables[ 5] = part2->Px();
858  fCandidateVariables[ 6] = part2->Py();
859  fCandidateVariables[ 7] = part2->Pz();
860  fCandidateVariables[ 8] = casc->MassXi();
861  fCandidateVariables[ 9] = TMath::Sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2));
862  if(xicobj->GetCharge()==-1) fCandidateVariables[10] = casc->MassLambda();
863  else fCandidateVariables[10] = casc->MassAntiLambda();
864  fCandidateVariables[11]= casc->Pt();
865  if(casc->ChargeXi()==-1) fCandidateVariables[12]= TMath::Sqrt(pow( casc->MomPosX(),2)+pow( casc->MomPosY(),2));
866  else fCandidateVariables[12]=TMath::Sqrt(pow( casc->MomNegX(),2)+pow( casc->MomNegY(),2));
867 
868  fCandidateVariables[13] = fVtx1->GetX();
869  fCandidateVariables[14] = fVtx1->GetY();
870  fCandidateVariables[15] = fVtx1->GetZ();
871  fCandidateVariables[16] = xicobj->GetOwnPrimaryVtx()->GetX();
872  fCandidateVariables[17] = xicobj->GetOwnPrimaryVtx()->GetY();
873  fCandidateVariables[18] = xicobj->GetOwnPrimaryVtx()->GetZ();
874 
875  fCandidateVariables[19] = xicobj->CascDcaXiDaughters();
876  fCandidateVariables[20] = xicobj->CascDcaV0Daughters();
877  fCandidateVariables[21] = xicobj->CascDecayLength();
883  fCandidateVariables[27] = xicobj->CascDecayLengthV0();
885 
886  fCandidateVariables[29] = dca[0];
887  fCandidateVariables[30] = dca[1];
888  fCandidateVariables[31] = dca[2];
889  fCandidateVariables[32] = xicobj->Getd0Prong(0);
890  fCandidateVariables[33] = xicobj->Getd0Prong(2);
891  fCandidateVariables[34] = xicobj->Getd0Prong(1);
892 
893  fCandidateVariables[35] = xicobj->DecayLength();
894  fCandidateVariables[36] = xicobj->DecayLengthXY();
897 
898  if(fAnalCuts->GetIsUsePID())
899  {
900  fCandidateVariables[39] = nSigmaTPCpi1;
901  fCandidateVariables[40] = nSigmaTPCpi2;
902  fCandidateVariables[41] = nSigmaTOFpi1;
903  fCandidateVariables[42] = nSigmaTOFpi2;
904  fCandidateVariables[43] = probPion1;
905  fCandidateVariables[44] = probPion2;
906  }
907 
908  fCandidateVariables[45] = -9999;
909  fCandidateVariables[46] = -9999;
910  fCandidateVariables[47] = -9999;
911  fCandidateVariables[48] = -9999;
912  fCandidateVariables[49] = -9999;
913  fCandidateVariables[50] = -9999;
914  fCandidateVariables[51] = -9999;
915  fCandidateVariables[52] = -9999;
916  fCandidateVariables[53] = -9999;
917  fCandidateVariables[54] = -9999;
918 
919 
920  if(fUseMCInfo){
921  if(mcpart){
922  fCandidateVariables[45] = mcpart->Label();
923  fCandidateVariables[46] = mcnused;
924  fCandidateVariables[47] = mcpart->GetPdgCode();
925  fCandidateVariables[51] = mcpart->Pt();
926  if(mcdaughter1&&mcdaughter2&&mcdaughterxi){
927  Double_t mcprimvertx = mcpart->Xv();
928  Double_t mcprimverty = mcpart->Yv();
929  Double_t mcsecvertx = mcdaughter1->Xv();
930  Double_t mcsecverty = mcdaughter1->Yv();
931  Double_t recosecvertx = xicobj->GetSecondaryVtx()->GetX();
932  Double_t recosecverty = xicobj->GetSecondaryVtx()->GetY();
933  fCandidateVariables[48] = TMath::Sqrt((mcsecvertx-mcprimvertx)*(mcsecvertx-mcprimvertx)+(mcsecverty-mcprimverty)*(mcsecverty-mcprimverty));
934  fCandidateVariables[49] = TMath::Sqrt((recosecvertx-mcprimvertx)*(recosecvertx-mcprimvertx)+(recosecverty-mcprimverty)*(recosecverty-mcprimverty));
935  Double_t vecx_vert = recosecvertx-mcprimvertx;
936  Double_t vecy_vert = recosecverty-mcprimverty;
937  Double_t vecl_vert = TMath::Sqrt(vecx_vert*vecx_vert+vecy_vert*vecy_vert);
938  Double_t vecx_mom = xicobj->Px();
939  Double_t vecy_mom = xicobj->Py();
940  Double_t vecl_mom = xicobj->Pt();
941  if(vecl_vert>0.&&vecl_mom>0.)
942  fCandidateVariables[50] = (vecx_vert*vecx_mom+vecy_vert*vecy_mom)/vecl_vert/vecl_mom;
943  fCandidateVariables[52] = mcdaughter1->Pt();
944  fCandidateVariables[53] = mcdaughter2->Pt();
945  fCandidateVariables[54] = mcdaughterxi->Pt();
946  }
947  }
948  }
949 
950  fCandidateVariables[55] = xicobj->GetSecondaryVtx()->GetX();
951  fCandidateVariables[56] = xicobj->GetSecondaryVtx()->GetY();
952  fCandidateVariables[57] = xicobj->GetSecondaryVtx()->GetZ();
953 
954  }//close if to check mc fill only signal
955  fVariablesTree->Fill();
956  }//fWriteTree
957 
958  if(fFillSparse){ //if you prefer to fill the histograms
959  Double_t trueXic=-9999;
960  if(isXic) trueXic=1.;
961  else trueXic=0.;
962  Double_t point[22]={xicobj->InvMassPiXiPi(), xicobj->Pt(), dca[2],dca[0], dca[1], xicobj->DecayLength(), xicobj->DecayLengthXY(), xicobj->XicCosPointingAngle(), xicobj->CascDcaXiDaughters(), xicobj->CascDcaV0Daughters(), xicobj->CascDcaV0ToPrimVertex(), xicobj->CascDcaPosToPrimVertex(), xicobj->CascDcaNegToPrimVertex(),xicobj->CascDcaBachToPrimVertex(),xicobj->CascCosPointingAngle(),xicobj->CascMassXi(), TMath::Sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)), nSigmaTPCpi1, nSigmaTPCpi2, nSigmaTOFpi1, nSigmaTOFpi2, trueXic};
963 
964  fSparseXicMass->Fill(point);
965  }
968  Double_t cont[3];
969  cont[0] = xicobj->InvMassPiXiPi();
970  cont[1] = xicobj->Pt();
971  cont[2] = fCentrality;
972  fHistoXicMass->Fill(cont);
973 
974  //these histos are filled only for candidates which pass IsSelected
975  fHistoDcaPi1Pi2->Fill(dca[2] );
976  fHistoDcaPi1Casc->Fill(dca[0]);
977  fHistoDcaPi2Casc->Fill(dca[1]);
978  fHistoLikeDecayLength->Fill(xicobj->DecayLength());
979  fHistoLikeDecayLengthXY->Fill(xicobj->DecayLengthXY());
980  fHistoXicCosPAXY->Fill(xicobj->XicCosPointingAngle());
981  fHistoXiMass->Fill(xicobj->CascMassXi());
989  fHistoXiPt->Fill(xicobj->PtProng(1));
990  fHistoPiPt->Fill(xicobj->PtProng(0));
991  fHistoPiPt->Fill(xicobj->PtProng(2));
992  fHistoPid0->Fill(xicobj->Getd0Prong(0));
993  fHistoPid0->Fill(xicobj->Getd0Prong(2));
994  fHistonSigmaTPCpi->Fill(nSigmaTPCpi1);
995  fHistonSigmaTPCpi->Fill(nSigmaTPCpi2);
996  fHistonSigmaTOFpi->Fill(nSigmaTOFpi1);
997  fHistonSigmaTOFpi->Fill(nSigmaTOFpi2);
998  fHistoProbPion->Fill(probPion1);
999  fHistoProbPion->Fill(probPion2);
1000  }
1001  }
1002  return;
1003 }
1004 
1005 //________________________________________________________________________
1007 {
1008  //
1009  // This is to define tree variables
1010  //
1011  const char* nameoutput = GetOutputSlot(3)->GetContainer()->GetName();
1012  fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
1013  Int_t nVar = 58;
1014  fCandidateVariables = new Float_t [nVar];
1015  TString * fCandidateVariableNames = new TString[nVar];
1016 
1017  fCandidateVariableNames[ 0]="InvMassXic";
1018  fCandidateVariableNames[ 1]="XicPt";
1019  fCandidateVariableNames[ 2]="Pi1Px";
1020  fCandidateVariableNames[ 3]="Pi1Py";
1021  fCandidateVariableNames[ 4]="Pi1Pz";
1022  fCandidateVariableNames[ 5]="Pi2Px";
1023  fCandidateVariableNames[ 6]="Pi2Py";
1024  fCandidateVariableNames[ 7]="Pi2Pz";
1025  fCandidateVariableNames[ 8]="MassXi"; //
1026  fCandidateVariableNames[ 9]="XiPt"; //
1027  fCandidateVariableNames[10]="MassLambda";
1028  fCandidateVariableNames[11]="LambdaPt";
1029  fCandidateVariableNames[12]="ProtonPt";
1030 
1031  fCandidateVariableNames[13]="PrimVtxX";
1032  fCandidateVariableNames[14]="PrimVtxY";
1033  fCandidateVariableNames[15]="PrimVtxZ";
1034  fCandidateVariableNames[16]="NewPrimVtxX";
1035  fCandidateVariableNames[17]="NewPrimVtxY";
1036  fCandidateVariableNames[18]="NewPrimVtxZ";
1037 
1038  fCandidateVariableNames[19]="CascDcaXiDaughters"; //
1039  fCandidateVariableNames[20]="CascDcaV0Daughters"; //
1040  fCandidateVariableNames[21]="CascDecayLengthXi";
1041  fCandidateVariableNames[22]="CascCosPointingAngleXi"; //
1042  fCandidateVariableNames[23]="CascDcaV0ToPrimVertex"; //
1043  fCandidateVariableNames[24]="CascDcaPosToPrimVertex"; //
1044  fCandidateVariableNames[25]="CascDcaNegToPrimVertex"; //
1045  fCandidateVariableNames[26]="CascDcaBachToPrimVertex"; //
1046  fCandidateVariableNames[27]="CascDecayLengthV0";
1047  fCandidateVariableNames[28]="CascCosPointingAngleV0";
1048 
1049  fCandidateVariableNames[29]="DcaPi1Casc"; //
1050  fCandidateVariableNames[30]="DcaPi2Casc"; //
1051  fCandidateVariableNames[31]="DcaPi1Pi2"; //
1052 
1053  fCandidateVariableNames[32]="Pi1d0"; //
1054  fCandidateVariableNames[33]="Pi2d0";
1055  fCandidateVariableNames[34]="Cascd0";
1056 
1057  fCandidateVariableNames[35]="DecayLength"; //for pions only?
1058  fCandidateVariableNames[36]="DecayLengthXY"; //for pions only??
1059  fCandidateVariableNames[37]="XicCosPAXY"; //
1060  fCandidateVariableNames[38]="BachelorsCosPAXY";
1061 
1062  fCandidateVariableNames[39]="nSigmaTPCpi1"; //
1063  fCandidateVariableNames[40]="nSigmaTPCpi2";
1064  fCandidateVariableNames[41]="nSigmaTOFpi1";//
1065  fCandidateVariableNames[42]="nSigmaTOFpi2";
1066  fCandidateVariableNames[43]="probPion1";
1067  fCandidateVariableNames[44]="probPion2";
1068 
1069  fCandidateVariableNames[45]="mcxicID";
1070  fCandidateVariableNames[46]="mcnused";
1071  fCandidateVariableNames[47]="mcpdgcode";
1072  fCandidateVariableNames[48]="mcdecaylength";
1073  fCandidateVariableNames[49]="mcdecaylength_secsmear";
1074  fCandidateVariableNames[50]="mcxiccospaxy";
1075  fCandidateVariableNames[51]="mcxicpt";
1076  fCandidateVariableNames[52]="mcpi1pt";
1077  fCandidateVariableNames[53]="mcpi2pt";
1078  fCandidateVariableNames[54]="mcxipt";
1079 
1080  fCandidateVariableNames[55]="XicVertX";
1081  fCandidateVariableNames[56]="XicVertY";
1082  fCandidateVariableNames[57]="XicVertZ";
1083 
1084  for (Int_t ivar=0; ivar<nVar; ivar++) {
1085  fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
1086  }
1087 
1088  return;
1089 }
1090 
1091 //__________________________________________________________________________
1093  //
1094  // This is to define general histograms
1095  //
1096 
1097  fCEvents = new TH1F("fCEvents","conter",25,-0.5,24.5);
1098  fCEvents->SetStats(kTRUE);
1099  fCEvents->GetXaxis()->SetBinLabel(1," ");
1100  fCEvents->GetXaxis()->SetBinLabel(2,"Analyzed events");
1101  fCEvents->GetXaxis()->SetBinLabel(3,"AliAODVertex exists");
1102  fCEvents->GetXaxis()->SetBinLabel(4,"Trigger not OK");
1103  fCEvents->GetXaxis()->SetBinLabel(5,"IsEventSelected");
1104  fCEvents->GetXaxis()->SetBinLabel(6,"CascadesHF exists");
1105  fCEvents->GetXaxis()->SetBinLabel(7,"MCarray exists");
1106  fCEvents->GetXaxis()->SetBinLabel(8,"MCheader exists");
1107  fCEvents->GetXaxis()->SetBinLabel(9,"triggerClass!=CINT1");
1108  fCEvents->GetXaxis()->SetBinLabel(10,"triggerMask!=kAnyINT");
1109  fCEvents->GetXaxis()->SetBinLabel(11,"triggerMask!=kAny");
1110  fCEvents->GetXaxis()->SetBinLabel(12,"vtxTitle.Contains(Z)");
1111  fCEvents->GetXaxis()->SetBinLabel(13,"vtxTitle.Contains(3D)");
1112  fCEvents->GetXaxis()->SetBinLabel(14,"vtxTitle.Doesn'tContain(Z-3D)");
1113  fCEvents->GetXaxis()->SetBinLabel(15,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
1114  fCEvents->GetXaxis()->SetBinLabel(16,"!IsEventSelected");
1115  fCEvents->GetXaxis()->SetBinLabel(17,"triggerMask!=kAnyINT || triggerClass!=CINT1");
1116  fCEvents->GetXaxis()->SetBinLabel(18,Form("zVtxMC<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
1117  fCEvents->GetXaxis()->SetBinLabel(19,"CascSelected");
1118  fCEvents->GetXaxis()->SetBinLabel(20,"#pi #pi pair selected");
1119  fCEvents->GetXaxis()->SetBinLabel(21,"correct charge combination");
1120  fCEvents->GetXaxis()->SetBinLabel(22,"#Xi_{c} with raugh sel");
1121  fCEvents->GetXaxis()->SetBinLabel(23,"Reco #pi #pi sec. vtx");
1122  fCEvents->GetXaxis()->SetBinLabel(24,"HF3Prongs");
1123 
1124  //fCEvents->GetXaxis()->SetTitle("");
1125  fCEvents->GetYaxis()->SetTitle("counts");
1126 
1127  fHTrigger = new TH1F("fHTrigger","counter",18,-0.5,17.5);
1128  fHTrigger->SetStats(kTRUE);
1129  fHTrigger->GetXaxis()->SetBinLabel(1," ");
1130  fHTrigger->GetXaxis()->SetBinLabel(2,"kMB");
1131  fHTrigger->GetXaxis()->SetBinLabel(3,"kSemiCentral");
1132  fHTrigger->GetXaxis()->SetBinLabel(4,"kCentral");
1133  fHTrigger->GetXaxis()->SetBinLabel(5,"kINT7");
1134  fHTrigger->GetXaxis()->SetBinLabel(6,"kEMC7");
1135  //fHTrigger->GetXaxis()->SetBinLabel(7,"Space");
1136  fHTrigger->GetXaxis()->SetBinLabel(8,"kMB|kSemiCentral|kCentral");
1137  fHTrigger->GetXaxis()->SetBinLabel(9,"kINT7|kEMC7");
1138  fHTrigger->GetXaxis()->SetBinLabel(11,"kMB&kSemiCentral");
1139  fHTrigger->GetXaxis()->SetBinLabel(12,"kMB&kCentral");
1140  fHTrigger->GetXaxis()->SetBinLabel(13,"kINT7&kEMC7");
1141 
1142  fHCentrality = new TH1F("fHCentrality","conter",100,0.,100.);
1143 
1144  Double_t binx[101];
1145  for(Int_t ib=0;ib<101;ib++){
1146  binx[ib] = 1.322-0.05 + 0.1/100.*(Double_t)ib ;
1147  }
1148  Double_t biny[21]={0.0,0.60,0.80,0.90,1.00,1.1,1.2,1.3,1.4,1.5,1.7,1.9,2.2,2.6,3.1,3.9,4.9,6.0,7.2,8.5,10.};
1149  fHistoXiMassvsPtRef1 = new TH2F("fHistoXiMassvsPtRef","Reference #Xi spectrum (-1.5 < #eta < -1.0)",100,binx,20,biny);
1150  fHistoXiMassvsPtRef2 = new TH2F("fHistoXiMassvsPtRef2","Reference #Xi spectrum (-1.0 < #eta < -0.5)",100,binx,20,biny);
1151  fHistoXiMassvsPtRef3 = new TH2F("fHistoXiMassvsPtRef3","Reference #Xi spectrum (-0.5 < #eta < -0.)",100,binx,20,biny);
1152  fHistoXiMassvsPtRef4 = new TH2F("fHistoXiMassvsPtRef4","Reference #Xi spectrum (0. < #eta < 0.5)",100,binx,20,biny);
1153  fHistoXiMassvsPtRef5 = new TH2F("fHistoXiMassvsPtRef5","Reference #Xi spectrum (0.5 < #eta < 1.)",100,binx,20,biny);
1154  fHistoXiMassvsPtRef6 = new TH2F("fHistoXiMassvsPtRef6","Reference #Xi spectrum (1. < #eta < 1.5)",100,binx,20,biny);
1155  fHistoPiPtRef = new TH1F("fHistoPiPtRef","Reference #pi spectrum",20,0.,10.);
1156  fHistoPiEtaRef = new TH1F("fHistoPiEtaRef","Reference #eta distributions of #pi ",50,-1,1.);
1157 
1158  fQAHistoNSelectedTracks = new TH1F("fQAHistoNSelectedTracks", "Number of tracks selected as pion candidates",100, 0, 100);
1159  fQAHistoNSelectedCasc = new TH1F("fQAHistoNSelectedCasc", "Number of tracks selected as cascades",20, 0, 20);
1160 
1161  fQAHistoDCApi1pi2 = new TH1F("fQAHistoDCApi1pi2","DCA #pi - #pi", 100,0.,0.5);
1162 
1163  fQAHistoAODPrimVertX = new TH1F("fQAHistoAODPrimVertX", "X coord of primary vertex", 100,0., 0.1);
1164  fQAHistoAODPrimVertY = new TH1F("fQAHistoAODPrimVertY", "Y coord of rimary vertex", 500,0., 0.5);
1165  fQAHistoAODPrimVertZ = new TH1F("fQAHistoAODPrimVertZ", "Z coord of primary vertex", 140,-12., 12.);
1166 
1167  fQAHistoRecoPrimVertX = new TH1F("fQAHistoRecoPrimVertX", "X coord of primary vertex without XiC decay tracks", 100,0., 0.1);
1168  fQAHistoRecoPrimVertY = new TH1F("fQAHistoRecoPrimVertY", "Y coord of primary vertex without XiC decay tracks", 500,0., 0.5);
1169  fQAHistoRecoPrimVertZ = new TH1F("fQAHistoRecoPrimVertZ", "Z coord of primary vertex without XiC decay tracks", 140,-12, 12);
1170 
1171  fQAHistoSecondaryVertexX = new TH1F("fQAHistoSecondaryVertexX", "X coord of secondary vertex", 1000, -2, 2);
1172  fQAHistoSecondaryVertexY = new TH1F("fQAHistoSecondaryVertexY", "Y coord of secondary vertex", 1000, -2, 2);
1173  fQAHistoSecondaryVertexZ = new TH1F("fQAHistoSecondaryVertexZ", "Z coord of secondary vertex", 1000, -20, 20);
1174  fQAHistoSecondaryVertexXY = new TH1F("fQAHistoSecondaryVertexXY", "XY coord of secondary vertex", 1000, -20, 20);
1175 
1176 
1177  fOutput->Add(fCEvents);
1178  fOutput->Add(fHTrigger);
1179  fOutput->Add(fHCentrality);
1186  fOutput->Add(fHistoPiPtRef);
1189  fOutput->Add(fQAHistoDCApi1pi2);
1200 
1201  if(fUseMCInfo) {
1202  fHistoMCSpectrumAccXic=new TH3F("fHistoMCSpectrumAccXic","fHistoMCSpectrumAccXic",250,0,50,20,-0.5,9.5,2,3.5,5.5);
1204  }
1205 
1206 
1207  return;
1208 }
1209 
1210 //__________________________________________________________________________
1212 {
1213  //
1214  // Define histograms
1215  //
1216 
1217  //------------------------------------------------
1218  // Basic histograms
1219  //------------------------------------------------
1220  Int_t bins_base[3]= {80,20,10};
1221  Double_t xmin_base[3]={2.468-0.2,0,0.00};
1222  Double_t xmax_base[3]={2.468+0.2,20.,100};
1223  fHistoXicMass = new THnSparseF("fHistoXicMass","fHistoXicMass;mass:pt:centrality",3,bins_base,xmin_base,xmax_base);
1224  fOutputAll->Add(fHistoXicMass);
1225 
1226  if (fFillSparse){
1227  Int_t bins_sparse[22]= { 100, 20, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 20, 50, 50, 50, 50, 50, 50, 2};
1228  Double_t xmin_sparse[22]={2.468-0.2, 0., 0., 0., 0., 0., 0., -1., 0., 0., 0., 0., 0., 0., 0.8, 1.322-0.05, 0., -5., -5., -5., -5., -0.5};
1229  Double_t xmax_sparse[22]={2.468+0.2, 20., 0.5, 0.5, 0.5, 0.2, 0.2, 1., 1., 1., 1., 1., 1., 1., 1., 1.322-0.05, 10., 5., 5., 5., 5., 1.5};
1230  fSparseXicMass = new THnSparseF("fSparseXicMass","fSparseXicMass;mass:pt:dcapi1pi2:dcapi1Xi:dcapi2Xi:decayLengthPi1Pi2:DecayLengthPi1P2XY:XicCosPAXY:DCAXidau:DCAV0dau:CascV0toPV:DCAPosToPrim:DCANegToPrim:DCABachToPrim:CosPAXiPrim:MassXi:XiPt:nSigmaTPCPi1:nSigmaTPCPi2:nSigmaTOFPi1:nSigmaTOFp2:isXic",22,bins_sparse,xmin_sparse, xmax_sparse);
1231  fOutputAll->Add(fSparseXicMass);
1232  }
1233  else {
1234  //------------------------------------------------
1235  // Checking histograms
1236  //------------------------------------------------
1237  fHistoDcaPi1Pi2 = new TH1F("fHistoDcaPi1Pi2","DCA (#pi_{1}-#pi_{2})",100,0.0,0.5);
1239  fHistoDcaPi1Casc = new TH1F("fHistoDcaPi1Casc","DCA (#pi_{1}-#Xi)",100,0.0,1.0);
1241  fHistoDcaPi2Casc = new TH1F("fHistoDcaPi2Casc","DCA (#pi_{2}-#Xi)",100,0.0,1.0);
1243  fHistoLikeDecayLength = new TH1F("fHistoLikeDecayLength","Decay Length (#pi-#pi)",100,0.,0.2);
1245  fHistoLikeDecayLengthXY = new TH1F("fHistoLikeDecayLengthXY","Decay Length (#pi-#pi)",100,0.,0.2);
1247  fHistoXicCosPAXY = new TH1F("fHistoXicCosPAXY","#Xi_{c} cos(pa) ",100,-1.0,1.0);
1249  fHistoXiMass=new TH1F("fHistoXiMass","#Xi^{-} Mass",100,1.322-0.05,1.322+0.05);
1250  fOutputAll->Add(fHistoXiMass);
1251  fHistoCascDcaXiDaughters=new TH1F("fHistoCascDcaXiDaughters","DCA #Xi daughters ",100,0.0,1.0);
1253  fHistoCascDcaV0Daughters=new TH1F("fHistoCascDcaV0Daughters","DCA #V0 daughters ",100,0.0,1.0);
1255  fHistoCascDcaV0ToPrimVertex=new TH1F("fHistoCascDcaV0ToPrimVertex","DCA V0 daughters to primary vertex ",100,0.0,1.0);
1257  fHistoCascDcaPosToPrimVertex=new TH1F("fHistoCascDcaPosToPrimVertex","DCA Pos-Prim ",100,0.0,1.0);
1259  fHistoCascDcaNegToPrimVertex=new TH1F("fHistoCascDcaNegToPrimVertex","DCA Neg-Prim ",100,0.0,1.0);
1261  fHistoCascDcaBachToPrimVertex=new TH1F("fHistoCascDcaBachToPrimVertex","DCA Bach-Prim ",100,0.0,1.0);
1263  fHistoCascCosPAXiPrim=new TH1F("fHistoCascCosPAXiPrim","#Xi CosPA (prim)",100,0.8,1.0);
1265  fHistoXiPt=new TH1F("fHistoXiPt","#Xi^{-} p_{T}",100,0.,10.);
1266  fOutputAll->Add(fHistoXiPt);
1267  fHistoPiPt=new TH1F("fHistoPiPt","#pi p_{T}",100,0.,10);
1268  fOutputAll->Add(fHistoPiPt);
1269  fHistoPid0=new TH1F("fHistoPid0","#pi d_{0}",100,-0.1,0.1);
1270  fOutputAll->Add(fHistoPid0);
1271  fHistonSigmaTPCpi=new TH1F("fHistonSigmaTPCpi","n#sigma (TPC, pion)",100,-10.,10.);
1273  fHistonSigmaTOFpi=new TH1F("fHistonSigmaTOFpi","n#sigma (TOF, pion)",100,-10.,10.);
1275  fHistoProbPion=new TH1F("fHistoProbPion","Bayse Prob",100,0.0,1.);
1276  fOutputAll->Add(fHistoProbPion);
1277  }
1278  return;
1279 }
1280 
1281 //________________________________________________________________________
1282 void AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::SelectTrack( const AliVEvent *event, Int_t trkEntries, Int_t &nSeleTrks,Bool_t *seleFlags)
1283 {
1284  //
1285  // Select good tracks using fAnalCuts (AliRDHFCuts object) and return the array of their ids
1286  //
1287 
1288  //const Int_t entries = event->GetNumberOfTracks();
1289  if(trkEntries==0) return;
1290 
1291  nSeleTrks=0;
1292  for(Int_t i=0; i<trkEntries; i++) {
1293  seleFlags[i] = kFALSE;
1294 
1295  AliVTrack *track;
1296  track = (AliVTrack*)event->GetTrack(i);
1297 
1298  if(track->GetID()<0) continue;
1299  Double_t covtest[21];
1300  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1301 
1302  AliAODTrack *aodt = (AliAODTrack*)track;
1303 
1304  if(!fAnalCuts) continue;
1305  if(fAnalCuts->SingleTrkCuts(aodt)){
1306  seleFlags[i]=kTRUE;
1307  nSeleTrks++;
1308  fHistoPiPtRef->Fill(aodt->Pt());
1309  fHistoPiEtaRef->Fill(aodt->Eta());
1310  }
1311  } // end loop on tracks
1312 }
1313 
1314 //________________________________________________________________________
1315 void AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::SelectCascade( const AliVEvent *event,Int_t nCascades,Int_t &nSeleCasc, Bool_t *seleCascFlags)
1316 {
1317  //
1318  // Select good cascade using fAnalCuts (AliRDHFCuts object) and return the array of their ids
1319  //
1320 
1321  Double_t primVtx[3];
1322  fVtx1->GetXYZ(primVtx);
1323 
1324  nSeleCasc = 0;
1325  for(Int_t icasc=0;icasc<nCascades;icasc++)
1326  {
1327  seleCascFlags[icasc] = kFALSE;
1328  AliAODcascade *casc = ((AliAODEvent*)event)->GetCascade(icasc);
1329 
1330  if(!fAnalCuts) continue;
1331  if(fAnalCuts->SingleCascadeCuts(casc,primVtx)){
1332  seleCascFlags[icasc] = kTRUE;
1333  nSeleCasc++;
1334  }
1335  if(fAnalCuts->SingleCascadeCutsRef(casc,primVtx)) //the only difference wrt SingleCascadeCuts is the cut on the Xi mass
1336  {
1337  Double_t rapxi = casc->RapXi();
1338  if(rapxi>-1.5&&rapxi<-1.0){
1339  fHistoXiMassvsPtRef1->Fill(casc->MassXi(),sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)));
1340  }
1341  if(rapxi>-1.0&&rapxi<-0.5){
1342  fHistoXiMassvsPtRef2->Fill(casc->MassXi(),sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)));
1343  }
1344  if(rapxi>-0.5&&rapxi<0.0){
1345  fHistoXiMassvsPtRef3->Fill(casc->MassXi(),sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)));
1346  }
1347  if(rapxi>0.0&&rapxi<0.5){
1348  fHistoXiMassvsPtRef4->Fill(casc->MassXi(),sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)));
1349  }
1350  if(rapxi>0.5&&rapxi<1.0){
1351  fHistoXiMassvsPtRef5->Fill(casc->MassXi(),sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)));
1352  }
1353  if(rapxi>1.0&&rapxi<1.5){
1354  fHistoXiMassvsPtRef6->Fill(casc->MassXi(),sqrt(pow(casc->MomXiX(),2)+pow(casc->MomXiY(),2)));
1355  }
1356  }
1357  } //loop over cascades
1358 }
1359 
1360 //________________________________________________________________________
1362 {
1363  //
1364  // Select LikeSign tracks
1365  //
1366 
1367  if(trk1->Charge()!=trk2->Charge()) return kFALSE;
1368  if(trk1->GetID()==trk2->GetID()) return kFALSE;
1369  return kTRUE;
1370 }
1371 
1372 //________________________________________________________________________
1373 AliAODVertex* AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::CallPrimaryVertex(AliAODcascade *casc, AliAODTrack *trk1, AliAODTrack *trk2, AliAODEvent* aod)
1374 {
1375  //
1376  // Make an array of tracks which should not be used in primary vertex calculation and
1377  // Call PrimaryVertex function
1378  //
1379 
1380  TObjArray *twoTrackArrayPlusXi = new TObjArray(5);
1381 
1382  AliESDtrack *cptrk1 = new AliESDtrack((AliVTrack*)trk1);
1383  twoTrackArrayPlusXi->AddAt(cptrk1,0);
1384  AliESDtrack *cptrk2 = new AliESDtrack((AliVTrack*)trk2);
1385  twoTrackArrayPlusXi->AddAt(cptrk2,1);
1386 
1387  AliESDtrack *cascptrack = new AliESDtrack((AliVTrack*)casc->GetDaughter(0));
1388  twoTrackArrayPlusXi->AddAt(cascptrack,2);
1389  AliESDtrack *cascntrack = new AliESDtrack((AliVTrack*)casc->GetDaughter(1));
1390  twoTrackArrayPlusXi->AddAt(cascntrack,3);
1391  AliESDtrack *cascbtrack = new AliESDtrack((AliVTrack*)casc->GetDecayVertexXi()->GetDaughter(0));
1392  twoTrackArrayPlusXi->AddAt(cascbtrack,4);
1393 
1394  AliAODVertex *newvert = PrimaryVertex(twoTrackArrayPlusXi,aod);
1395 
1396  for(Int_t i=0;i<5;i++)
1397  {
1398  AliESDtrack *tesd = (AliESDtrack*)twoTrackArrayPlusXi->UncheckedAt(i);
1399  delete tesd;
1400  }
1401  twoTrackArrayPlusXi->Clear();
1402  delete twoTrackArrayPlusXi;
1403 
1404  return newvert;
1405 }
1406 
1407 //________________________________________________________________________
1408 AliAODVertex* AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::CallReconstructSecondaryVertex(AliAODTrack *trk1, AliAODTrack *trk2, Double_t &dispersion)
1409 {
1410  //
1411  // Make an array of tracks which I want to use for vertex calculation and call ReconstructSecondaryVertex
1412  //
1413 
1414  TObjArray *trkArray;
1415  trkArray= new TObjArray(2);
1416 
1417  AliESDtrack *cptrk1 = new AliESDtrack((AliVTrack*)trk1);
1418  trkArray->AddAt(cptrk1,0);
1419  AliESDtrack *cptrk2 = new AliESDtrack((AliVTrack*)trk2);
1420  trkArray->AddAt(cptrk2,1);
1421 
1422  Double_t xdummy, ydummy;
1423  Double_t dcap1p2 = cptrk1->GetDCA(cptrk2,fBzkG,xdummy,ydummy);
1424  fQAHistoDCApi1pi2->Fill(dcap1p2);
1425  AliAODVertex *secvertex=0;
1426  if(dcap1p2<fAnalCuts->GetProdLikeSignDcaMax()){
1427  secvertex = ReconstructSecondaryVertex(trkArray,dispersion);
1428  }
1429 
1430  for(Int_t i=0;i<2;i++)
1431  {
1432  AliESDtrack *tesd = (AliESDtrack*)trkArray->UncheckedAt(i);
1433  delete tesd;
1434  }
1435  trkArray->Clear();
1436  delete trkArray;
1437 
1438  return secvertex;
1439 }
1440 
1441 //________________________________________________________________________
1443  AliVEvent *event)
1444 {
1445  //
1446  //Used only for pp
1447  //copied from AliAnalysisVertexingHF (except for the following 3 lines)
1448  //
1449 
1450  Bool_t fRecoPrimVtxSkippingTrks = kTRUE;
1451  Bool_t fRmTrksFromPrimVtx = kFALSE;
1452 
1453  AliESDVertex *vertexESD = 0;
1454  AliAODVertex *vertexAOD = 0;
1455 
1456  //vertexESD = new AliESDVertex(*fV1);
1457 
1458 
1459  if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1460  // primary vertex from the input event
1461 
1462  vertexESD = new AliESDVertex(*fV1);
1463 
1464  } else {
1465  // primary vertex specific to this candidate
1466 
1467  Int_t nTrks = trkArray->GetEntriesFast();
1468  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1469 
1470  if(fRecoPrimVtxSkippingTrks) {
1471  // recalculating the vertex
1472 
1473  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1474  Float_t diamondcovxy[3];
1475  event->GetDiamondCovXY(diamondcovxy);
1476  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1477  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1478  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1479  vertexer->SetVtxStart(diamond);
1480  delete diamond; diamond=NULL;
1481  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1482  vertexer->SetOnlyFitter();
1483  }
1484  Int_t skipped[1000];
1485  Int_t nTrksToSkip=0,id;
1486  AliExternalTrackParam *t = 0;
1487  for(Int_t i=0; i<nTrks; i++) {
1488  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1489  id = (Int_t)t->GetID();
1490  if(id<0) continue;
1491  skipped[nTrksToSkip++] = id;
1492  }
1493  // TEMPORARY FIX
1494  // For AOD, skip also tracks without covariance matrix
1495  Double_t covtest[21];
1496  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1497  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1498  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1499  id = (Int_t)vtrack->GetID();
1500  if(id<0) continue;
1501  skipped[nTrksToSkip++] = id;
1502  }
1503  }
1504  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
1505  //
1506  vertexer->SetSkipTracks(nTrksToSkip,skipped);
1507  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1508 
1509  } else if(fRmTrksFromPrimVtx && nTrks>0) {
1510  // removing the prongs tracks
1511 
1512  TObjArray rmArray(nTrks);
1513  UShort_t *rmId = new UShort_t[nTrks];
1514  AliESDtrack *esdTrack = 0;
1515  AliESDtrack *t = 0;
1516  for(Int_t i=0; i<nTrks; i++) {
1517  t = (AliESDtrack*)trkArray->UncheckedAt(i);
1518  esdTrack = new AliESDtrack(*t);
1519  rmArray.AddLast(esdTrack);
1520  if(esdTrack->GetID()>=0) {
1521  rmId[i]=(UShort_t)esdTrack->GetID();
1522  } else {
1523  rmId[i]=9999;
1524  }
1525  }
1526  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
1527  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1528  delete [] rmId; rmId=NULL;
1529  rmArray.Delete();
1530 
1531  }
1532 
1533  delete vertexer; vertexer=NULL;
1534  if(!vertexESD) return vertexAOD;
1535  if(vertexESD->GetNContributors()<=0) {
1536  //AliDebug(2,"vertexing failed");
1537  delete vertexESD; vertexESD=NULL;
1538  return vertexAOD;
1539  }
1540 
1541 
1542  }
1543 
1544  // convert to AliAODVertex
1545  Double_t pos[3],cov[6],chi2perNDF;
1546  vertexESD->GetXYZ(pos); // position
1547  vertexESD->GetCovMatrix(cov); //covariance matrix
1548  chi2perNDF = vertexESD->GetChi2toNDF();
1549  delete vertexESD; vertexESD=NULL;
1550 
1551  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1552 
1553  return vertexAOD;
1554 }
1555 
1556 //________________________________________________________________________
1558 {
1559  //
1560  // Reconstruct secondary vertex from trkArray (Copied from AliAnalysisVertexingHF)
1561  //
1562 
1563  AliESDVertex *vertexESD = 0;
1564  AliAODVertex *vertexAOD = 0;
1565 
1566  AliVertexerTracks *fVertexerTracks = new AliVertexerTracks(fBzkG);
1567 
1568  fVertexerTracks->SetVtxStart(fV1);
1569  vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
1570 
1571  delete fVertexerTracks; fVertexerTracks=NULL;
1572 
1573  if(!vertexESD) return vertexAOD;
1574 
1575  if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
1576  //AliDebug(2,"vertexing failed");
1577  delete vertexESD; vertexESD=NULL;
1578  return vertexAOD;
1579  }
1580 
1581  Double_t vertRadius2=vertexESD->GetX()*vertexESD->GetX()+vertexESD->GetY()*vertexESD->GetY();
1582  if(vertRadius2>8.){
1583  // vertex outside beam pipe, reject candidate to avoid propagation through material
1584  delete vertexESD; vertexESD=NULL;
1585  return vertexAOD;
1586  }
1587 
1588 
1589  // convert to AliAODVertex
1590  Double_t pos[3],cov[6],chi2perNDF;
1591  vertexESD->GetXYZ(pos); // position
1592  vertexESD->GetCovMatrix(cov); //covariance matrix
1593  chi2perNDF = vertexESD->GetChi2toNDF();
1594  dispersion = vertexESD->GetDispersion();
1595  delete vertexESD; vertexESD=NULL;
1596 
1597  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
1598  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1599 
1600  return vertexAOD;
1601 }
1602 
1603 //________________________________________________________________________
1604 AliAODRecoCascadeHF3Prong* AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::MakeCascadeHF3Prong(AliAODcascade *casc, AliAODTrack *part1, AliAODTrack *part2, AliAODEvent * aod, AliAODVertex *secVert, Double_t dispersion)
1605 {
1606  //
1607  // Make AliAODRecoCascadeHF3Prong object from the arguments
1608  //
1609 
1610  //------------------------------------------------
1611  // PrimaryVertex
1612  //------------------------------------------------
1613  AliAODVertex *primVertexAOD;
1614  Bool_t unsetvtx = kFALSE;
1616  primVertexAOD = CallPrimaryVertex(casc,part1,part2,aod);
1617  if(!primVertexAOD){
1618  primVertexAOD = fVtx1;
1619  }else{
1620  unsetvtx = kTRUE;
1621  }
1622  }else{
1623  primVertexAOD = fVtx1;
1624  }
1625  if(!primVertexAOD) return 0x0;
1626  Double_t pos[3]; primVertexAOD->GetXYZ(pos);
1627  fQAHistoRecoPrimVertX->Fill(pos[0]);
1628  fQAHistoRecoPrimVertY->Fill(pos[1]);
1629  fQAHistoRecoPrimVertZ->Fill(pos[2]);
1630 
1631 
1632  //------------------------------------------------
1633  // DCA between tracks
1634  //------------------------------------------------
1635  AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)part1);
1636  AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)part2);
1637  Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
1638  xyz[0]=casc->DecayVertexXiX();
1639  xyz[1]=casc->DecayVertexXiY();
1640  xyz[2]=casc->DecayVertexXiZ();
1641  pxpypz[0]=casc->MomXiX();
1642  pxpypz[1]=casc->MomXiY();
1643  pxpypz[2]=casc->MomXiZ();
1644  casc->GetCovarianceXYZPxPyPz(cv);
1645  sign=casc->ChargeXi();
1646  AliExternalTrackParam *trackCasc = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
1647 
1648  Double_t xdummy, ydummy;
1649  Double_t dcap1p2 = esdtrack1->GetDCA(esdtrack2,fBzkG,xdummy,ydummy);
1650  Double_t dcap1casc = esdtrack1->GetDCA(trackCasc,fBzkG,xdummy,ydummy);
1651  Double_t dcap2casc = esdtrack2->GetDCA(trackCasc,fBzkG,xdummy,ydummy);
1652  Double_t dca[3]={dcap1casc,dcap2casc,dcap1p2};
1653 
1654  //------------------------------------------------
1655  // Propagate all tracks to the secondary vertex and calculate momentum there
1656  //------------------------------------------------
1657 
1658  Double_t d0z0[2],covd0z0[3];
1659  Double_t secR = TMath::Sqrt(secVert->GetX()*secVert->GetX()+secVert->GetY()*secVert->GetY());
1660  if(secR<1.0){
1661  part1->PropagateToDCA(secVert,fBzkG,kVeryBig,d0z0,covd0z0);
1662  part2->PropagateToDCA(secVert,fBzkG,kVeryBig,d0z0,covd0z0);
1663  trackCasc->PropagateToDCA(secVert,fBzkG,kVeryBig);
1664  }else{
1665  part1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1666  part2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1667  trackCasc->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig);
1668  }
1669 
1670  Double_t momxi_new[3]={-9999,-9999,-9999.};
1671  trackCasc->GetPxPyPz(momxi_new);
1672 
1673  Double_t px[3],py[3],pz[3];
1674  px[0] = part1->Px(); py[0] = part1->Py(); pz[0] = part1->Pz();
1675  px[1] = momxi_new[0]; py[1] = momxi_new[1]; pz[1] = momxi_new[2];
1676  px[2] = part2->Px(); py[2] = part2->Py(); pz[2] = part2->Pz();
1677 
1678  //------------------------------------------------
1679  // d0
1680  //------------------------------------------------
1681  Double_t d0[3],d0err[3];
1682  part1->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1683  d0[0]= d0z0[0];
1684  d0err[0] = TMath::Sqrt(covd0z0[0]);
1685 
1686  trackCasc->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1687  d0[1]= d0z0[0];
1688  d0err[1] = TMath::Sqrt(covd0z0[0]);//Do not know what to use
1689 
1690  part2->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1691  d0[2]= d0z0[0];
1692  d0err[2] = TMath::Sqrt(covd0z0[0]);
1693 
1694 
1695  //------------------------------------------------
1696  // Other stuff
1697  //------------------------------------------------
1698  Double_t dist12=0.0;//Vertex pi-casc not calculated
1699  Double_t dist23=0.0;//Vertex pi-casc not calculated
1700  Short_t charge=(Short_t)(part1->Charge()+trackCasc->Charge()+part2->Charge());
1701 
1702  //------------------------------------------------
1703  // Create AliAODRecoCascadeHF3Prong
1704  //------------------------------------------------
1705  AliAODRecoCascadeHF3Prong *theCascade = new AliAODRecoCascadeHF3Prong(secVert,charge,px,py,pz,d0,d0err,dca,dispersion,dist12,dist23);
1706  if(!theCascade)
1707  {
1708  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1709  if(esdtrack1) delete esdtrack1;
1710  if(esdtrack2) delete esdtrack2;
1711  if(trackCasc) delete trackCasc;
1712  return 0x0;
1713  }
1714  theCascade->SetOwnPrimaryVtx(primVertexAOD);
1715  UShort_t id[3]={(UShort_t)part1->GetID(),(UShort_t)trackCasc->GetID(),(UShort_t)part2->GetID()};
1716  theCascade->SetProngIDs(3,id);
1717 
1718  theCascade->GetSecondaryVtx()->AddDaughter(part1);
1719  theCascade->GetSecondaryVtx()->AddDaughter(casc);
1720  theCascade->GetSecondaryVtx()->AddDaughter(part2);
1721 
1722  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1723  if(esdtrack1) delete esdtrack1;
1724  if(esdtrack2) delete esdtrack2;
1725  if(trackCasc) delete trackCasc;
1726 
1727  return theCascade;
1728 }
1729 
1730 //________________________________________________________________________
1732  cout<< "Executing Loop over Gen Particles\n";
1733  for(Int_t kmc=0; kmc<mcArray->GetEntries(); kmc++){
1734  AliAODMCParticle *mcpart=(AliAODMCParticle*)mcArray->At(kmc);
1735  if(!mcpart) continue;
1736 
1737  Int_t pdg=TMath::Abs(mcpart->GetPdgCode());
1738  Int_t arrayDauLab[5];
1739  if(pdg==4232){
1740  if(CheckXic2XiPiPi(mcArray,mcpart,arrayDauLab)==1){ //the arrayDauLab is used to check if the single particles are in acceptance when the Xic is in the acceptance.
1741  Int_t checkOrigin=AliVertexingHFUtils::CheckOrigin(mcArray,mcpart,kTRUE);
1742  if(checkOrigin==0)continue;
1743 
1744  Float_t ptpart=mcpart->Pt();
1745  Float_t ypart=mcpart->Y();
1746 
1747  if(TMath::Abs(ypart)<0.5){
1748  fHistoMCSpectrumAccXic->Fill(ptpart,kGenLimAcc,checkOrigin);
1749  }
1750  Bool_t isInAcc=kTRUE;
1751  // check GenAcc level
1752  if(fAnalCuts){
1753  if(!fAnalCuts->IsInFiducialAcceptance(ptpart,ypart)){
1754  isInAcc=kFALSE;
1755  }
1756  } else if (TMath::Abs(ypart)>0.8) isInAcc=kFALSE;
1757 
1758  if(isInAcc){
1759  fHistoMCSpectrumAccXic->Fill(ptpart,kGenAccMother,checkOrigin);
1760  for(Int_t k=0;k<5;k++){
1761  AliAODMCParticle *mcpartdau=(AliAODMCParticle*)mcArray->At(arrayDauLab[k]);
1762  if(TMath::Abs(mcpartdau->Eta())>0.9) isInAcc=kFALSE;
1763  }
1764  }
1765  if(isInAcc) fHistoMCSpectrumAccXic->Fill(ptpart,kGenAcc,checkOrigin);
1766  } //else continue;//CheckXic2XiPiPi
1767  }//Check on PDG code
1768  }//loop over particles
1769 }
1770 
1771 //________________________________________________________________________
1772 Int_t AliAnalysisTaskSEXicPlus2XiPiPifromAODtracks::CheckXic2XiPiPi(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1773  cout<< "Executing Xic2XiPiPi\n";
1774  if(!mcPart) return -1;
1775 
1776  Int_t pdgD=mcPart->GetPdgCode();
1777 
1778  if(TMath::Abs(pdgD)!=4232) return -1;
1779 
1780  Int_t nDau=mcPart->GetNDaughters();
1781  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1782 
1783  Int_t nXi=0;
1784  Int_t nLambda=0;
1785  Int_t nPions=0;
1786  Int_t nProtons=0;
1787  Double_t sumPxDau=0.;
1788  Double_t sumPyDau=0.;
1789  Double_t sumPzDau=0.;
1790  Int_t nFound=0;
1791 
1792  Int_t codeRes=-1;
1793  if(nDau==3){
1794  for(Int_t iDau=0; iDau<nDau; iDau++){
1795  cout << "XiC daughter " << iDau << endl;
1796  Int_t indDau = labelFirstDau+iDau;
1797  if(indDau<0) return -1;
1798  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1799  if(!dau) return -1;
1800  Int_t pdgdau=dau->GetPdgCode();
1801  if(TMath::Abs(pdgdau)==3312){ //Xi
1802  cout << "I am a Xi" << endl;
1803  nXi++;
1804  sumPxDau+=dau->Px();
1805  sumPyDau+=dau->Py();
1806  sumPzDau+=dau->Pz();
1807  Int_t nXiDau=dau->GetNDaughters();
1808  if(nXiDau!=2) return -1; //no correct Xi decay
1809  Int_t indFirstXiDau=dau->GetDaughterLabel(0);
1810  for(Int_t XiDau=0; XiDau<2; XiDau++){
1811  Int_t indXiDau=indFirstXiDau+XiDau;
1812  if(indXiDau<0) return -1;
1813  AliAODMCParticle* Xidau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indXiDau));
1814  if(!Xidau) return -1;
1815  Int_t pdgXidau=Xidau->GetPdgCode();
1816  if(TMath::Abs(pdgXidau)==3122){
1817  cout<< "-- I am a lambda" <<endl;
1818  nLambda++;
1819  Int_t nLambdaDau=Xidau->GetNDaughters();
1820  if (nLambdaDau!=2) return -1;
1821  Int_t indFirstLambdaDau=Xidau->GetDaughterLabel(0);
1822  for(Int_t LambdaDau=0; LambdaDau<2; LambdaDau++){
1823  Int_t indLambdaDau=indFirstLambdaDau+LambdaDau;
1824  AliAODMCParticle* Lambdadau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indLambdaDau));
1825  if(!Lambdadau) return -1;
1826  Int_t pdgLambdadau=Lambdadau->GetPdgCode();
1827  if(TMath::Abs(pdgLambdadau)==2212){
1828  cout<< "-- I am a proton from Lambda" <<endl;
1829  nProtons++;
1830  cout << nFound << endl;
1831  arrayDauLab[nFound++]=indLambdaDau;
1832  if (nFound>5) return -1;
1833  } else if(TMath::Abs(pdgLambdadau)==211){
1834  cout<< "-- I am a pion from Lambda" <<endl;
1835  nPions++;
1836  cout << nFound << endl;
1837  arrayDauLab[nFound++]=indLambdaDau;
1838  if (nFound>5) return -1;
1839  }
1840  }
1841  }else if(TMath::Abs(pdgXidau)==211) {
1842  cout<< "-- I am a pion from Xi" <<endl;
1843  nPions++;
1844  cout << nFound << endl;
1845  arrayDauLab[nFound++]=indXiDau;
1846  if (nFound>5) return -1;
1847  }
1848  }
1849  }else if(TMath::Abs(pdgdau)==211){
1850  cout << "I am a pion from Xic " << endl;
1851  nPions++;
1852  sumPxDau+=dau->Px();
1853  sumPyDau+=dau->Py();
1854  sumPzDau+=dau->Pz();
1855  cout << nFound << endl;
1856  arrayDauLab[nFound++]=indDau;
1857  if (nFound>5) return -1;
1858  }
1859  }
1860  if(nXi!=1) return -1;
1861  if(nLambda!=1) return -1;
1862  if(nPions!=4) return -1;
1863  if(nProtons!=1) return -1;
1864  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1865  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1866  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1867 
1868  if(nXi==1 && nLambda==1 && nPions==4 && nProtons ==1) return 1;
1869  } else return -1;
1870 }
1871 
Int_t charge
Double_t InvMassPiXiPi() const
Xic invariant mass.
Int_t pdg
AliAODTrack * GetBachelor2() const
double Double_t
Definition: External.C:58
Definition: External.C:260
Bool_t SingleCascadeCutsRef(AliAODcascade *casc, Double_t *vert)
AliAODVertex * PrimaryVertex(const TObjArray *trkArray, AliVEvent *event)
void SelectTrack(const AliVEvent *event, Int_t trkEntries, Int_t &nSeleTrks, Bool_t *seleFlags)
Definition: External.C:236
Bool_t GetUseCombined()
Definition: AliAODPidHF.h:178
char Char_t
Definition: External.C:18
TH1F * fQAHistoNSelectedCasc
! QA histo for number of selected Cascades/event
AliAODVertex * CallReconstructSecondaryVertex(AliAODTrack *trk1, AliAODTrack *trk2, Double_t &disp)
TList * fOutputAll
! User output slot 3 // Analysis histos
TH1F * fQAHistoSecondaryVertexX
! Coordinates of the reconstructed secondary vertex
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
void SelectCascade(const AliVEvent *event, Int_t nCascades, Int_t &nSeleCasc, Bool_t *seleCascFlags)
TH1F * fQAHistoRecoPrimVertZ
! Coordinates of the reconstructed primary vertex without XiC decay tracks
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
AliAODcascade * GetCascade() const
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:274
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:264
TH1F * fQAHistoSecondaryVertexZ
! Coordinates of the reconstructed secondary vertex
Bool_t fIsMB
Reconstruct primary vertex excluding candidate tracks.
AliAODVertex * ReconstructSecondaryVertex(TObjArray *trkArray, Double_t &dispersion, Bool_t useTRefArray=kTRUE)
THnSparse * fSparseXicMass
! xic sparse to study cut variation
AliNormalizationCounter * fCounter
!Counter for normalization slot 4
TH1F * fQAHistoRecoPrimVertX
! Coordinates of the reconstructed primary vertex without XiC decay tracks
virtual void UserCreateOutputObjects()
Implementation of interface methods.
int Int_t
Definition: External.C:63
Int_t CheckXic2XiPiPi(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Int_t *arrayDauLab)
float Float_t
Definition: External.C:68
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t fFillSparse
flag to decide whether to write the candidate variables on a tree variables
void FillROOTObjects(AliAODRecoCascadeHF3Prong *xicobj, AliAODMCParticle *mcpart, AliAODMCParticle *mcdau1, AliAODMCParticle *mcdau2, AliAODMCParticle *mcdauxi, Int_t mcnused, Bool_t isXiC)
AliAODTrack * GetBachelor1() const
AliAODVertex * GetOwnPrimaryVtx() const
TH1F * fQAHistoRecoPrimVertY
! Coordinates of the reconstructed primary vertex without XiC decay tracks
AliAODVertex * CallPrimaryVertex(AliAODcascade *casc, AliAODTrack *trk1, AliAODTrack *trk2, AliAODEvent *evt)
void SetProngIDs(Int_t nIDs, UShort_t *id)
TH1F * fHistoCascDcaPosToPrimVertex
! DCA of positive track to primary vertex
AliPIDResponse * GetPidResponse() const
Definition: AliAODPidHF.h:173
void GetDCAs(Double_t dca[3]) const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
short Short_t
Definition: External.C:23
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Double_t DecayLengthXY() const
TH1F * fHistoCascDcaBachToPrimVertex
! DCA of bachelor track to primary vertex
AliAODRecoCascadeHF3Prong * MakeCascadeHF3Prong(AliAODcascade *casc, AliAODTrack *trk1, AliAODTrack *trk2, AliAODEvent *aod, AliAODVertex *secvert, Double_t dispersion)
Bool_t SingleCascadeCuts(AliAODcascade *casc, Double_t *vert)
TH1F * fQAHistoNSelectedTracks
! QA histo for number of selected tracks/event
TH1F * fQAHistoSecondaryVertexY
! Coordinates of the reconstructed secondary vertex
AliRDHFCutsXicPlustoXiPiPifromAODtracks * fAnalCuts
histogram to check centrality
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:221
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabscasc, Int_t *pdgDg, Int_t *pdgDgcasc, Int_t *pdgDgv0, TClonesArray *mcArray) const
TH1F * fQAHistoSecondaryVertexXY
! Coordinates of the reconstructed secondary vertex
TH1F * fHistoCascCosPAXiPrim
! Cosine pointing angle of Xi to primary vertex
unsigned short UShort_t
Definition: External.C:28
Bool_t SelectWithRoughCuts(AliAODcascade *casc, AliAODTrack *trk1, AliAODTrack *trk2)
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:287
const char Option_t
Definition: External.C:48
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:343
bool Bool_t
Definition: External.C:53
TH1F * fHistoCascDcaNegToPrimVertex
! DCA of negative track to primary vertex
Double_t DecayLength() const
TH1F * fQAHistoDCApi1pi2
! QA histo for dca betwen two pions from XiC
TTree * fVariablesTree
flag to decide whether fill the THnSparse