AliPhysics  ced2227 (ced2227)
AliAnalysisTaskSELc2pK0sfromAODtracks.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 // Lc->pK0s 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 // Use "AddCutStream" function to add a cut.
30 //
31 //-------------------------------------------------------------------------
32 //
33 // Authors: Y.S Watanabe(a)
34 // (a) CNS, the University of Tokyo
35 // Contatcs: wyosuke@cns.s.u-tokyo.ac.jp
36 //-------------------------------------------------------------------------
37 
38 #include <TSystem.h>
39 #include <TParticle.h>
40 #include <TParticlePDG.h>
41 #include <TH1F.h>
42 #include <TH1F.h>
43 #include <TH2F.h>
44 #include <TH3F.h>
45 #include <TRandom.h>
46 #include <THnSparse.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"
63 #include "AliAODRecoCascadeHF.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 "AliTOFPIDResponse.h"
74 #include "AliAODPidHF.h"
75 #include "AliInputEventHandler.h"
76 #include "AliESDtrackCuts.h"
77 #include "AliNeutralTrackParam.h"
78 #include "AliKFParticle.h"
79 #include "AliKFVertex.h"
80 #include "AliExternalTrackParam.h"
81 #include "AliESDtrack.h"
82 #include "AliCentrality.h"
83 #include "AliVertexerTracks.h"
85 #include <set>
86 
87 using std::cout;
88 using std::endl;
89 
93 
94 //__________________________________________________________________________
97  fUseMCInfo(kFALSE),
98  fOutput(0),
99  fOutputAll(0),
100  fListCuts(0),
101  fCEvents(0),
102  fHTrigger(0),
103  fHCentrality(0),
104  fHReactionPlane(0),
105  fAnalCuts(0),
106  fIsEventSelected(kFALSE),
107  fWriteVariableTree(kFALSE),
108  fWriteEachVariableTree(kFALSE),
109  fWriteMCVariableTree(kFALSE),
110  fVariablesTree(0),
111  fProtonVariablesTree(0),
112  fV0VariablesTree(0),
113  fMCVariablesTree(0),
114  fMCProtonVariablesTree(0),
115  fMCV0VariablesTree(0),
116  fReconstructPrimVert(kFALSE),
117  fIsMB(kFALSE),
118  fIsSemi(kFALSE),
119  fIsCent(kFALSE),
120  fIsINT7(kFALSE),
121  fIsEMC7(kFALSE),
122  fCandidateVariables(),
123  fCandidateProtonVariables(),
124  fCandidateV0Variables(),
125  fCandidateMCVariables(),
126  fCandidateMCProtonVariables(),
127  fCandidateMCV0Variables(),
128  fVtx1(0),
129  fV1(0),
130  fVtxZ(0),
131  fBzkG(0),
132  fCentrality(0),
133  fReactionPlane(0),
134  fRunNumber(0),
135  fTriggerCheck(0),
136  fUseCentralityV0M(kFALSE),
137  fEvNumberCounter(0),
138  fCounter(0),
139  fHistonEvtvsRunNumber(0),
140  fHistonProtonvsRunNumber(0),
141  fHistonK0svsRunNumber(0),
142  fHistoLcMCGen(0),
143  fHistoLcK0SpMass(0),
144  fHistoLcK0SpMassMix(0),
145  fHistoLcK0SpMassCoarse(0),
146  fHistoLcK0SpMassMixCoarse(0),
147  fHistoK0spCorrelation(0),
148  fHistoK0spCorrelationMix(0),
149  fHistoK0spCorrelationMCS(0),
150  fHistoLcK0SpMassMCS(0),
151  fHistoLcK0SpPi0MassMCS(0),
152  fHistoLcKPluspMass(0),
153  fHistoLcKMinuspMass(0),
154  fHistoLcKPluspMassMix(0),
155  fHistoLcKMinuspMassMix(0),
156  fHistoBachPt(0),
157  fHistoBachPtMCS(0),
158  fHistoBachPtMCGen(0),
159  fHistoKaonPt(0),
160  fHistoKaonPtMCS(0),
161  fHistoKaonPtMCGen(0),
162  fHistoK0sMassvsPt(0),
163  fHistoK0sMassvsPtMCS(0),
164  fHistoK0sMassvsPtMCGen(0),
165  fHistod0Bach(0),
166  fHistod0V0(0),
167  fHistod0d0(0),
168  fHistoV0CosPA(0),
169  fHistoProbProton(0),
170  fHistoDecayLength(0),
171  fHistoK0SMass(0),
172  fHistoMassTagV0Min(0),
173  fHistoMassTagV0SameSignMin(0),
174  fHistoResponseLcPt(0),
175  fHistoResponseLcPt1(0),
176  fHistoResponseLcPt2(0),
177  fGTI(0),fGTIndex(0), fTrackBuffSize(19000),
178  fDoEventMixing(0),
179  fNumberOfEventsForMixing (5),
180  fNzVtxBins (0),
181  fNCentBins (0),
182  fNRPBins (0),
183  fNOfPools(1),
184  fEventBuffer(0x0),
185  fEventInfo(0x0),
186  fProtonTracks(0x0),
187  fV0Tracks(0x0),
188  fProtonCutVarsArray(0x0),
189  fV0CutVarsArray(0x0)
190 {
191  //
192  // Default Constructor.
193  //
194 }
195 
196 //___________________________________________________________________________
199  Bool_t writeVariableTree) :
200  AliAnalysisTaskSE(name),
201  fUseMCInfo(kFALSE),
202  fOutput(0),
203  fOutputAll(0),
204  fListCuts(0),
205  fCEvents(0),
206  fHTrigger(0),
207  fHCentrality(0),
208  fHReactionPlane(0),
209  fAnalCuts(analCuts),
210  fIsEventSelected(kFALSE),
211  fWriteVariableTree(writeVariableTree),
212  fWriteEachVariableTree(kFALSE),
213  fWriteMCVariableTree(kFALSE),
214  fVariablesTree(0),
216  fV0VariablesTree(0),
217  fMCVariablesTree(0),
220  fReconstructPrimVert(kFALSE),
221  fIsMB(kFALSE),
222  fIsSemi(kFALSE),
223  fIsCent(kFALSE),
224  fIsINT7(kFALSE),
225  fIsEMC7(kFALSE),
232  fVtx1(0),
233  fV1(0),
234  fVtxZ(0),
235  fBzkG(0),
236  fCentrality(0),
237  fReactionPlane(0),
238  fRunNumber(0),
239  fTriggerCheck(0),
240  fUseCentralityV0M(kFALSE),
241  fEvNumberCounter(0),
242  fCounter(0),
246  fHistoLcMCGen(0),
247  fHistoLcK0SpMass(0),
260  fHistoBachPt(0),
261  fHistoBachPtMCS(0),
263  fHistoKaonPt(0),
264  fHistoKaonPtMCS(0),
269  fHistod0Bach(0),
270  fHistod0V0(0),
271  fHistod0d0(0),
272  fHistoV0CosPA(0),
273  fHistoProbProton(0),
275  fHistoK0SMass(0),
281  fGTI(0),fGTIndex(0), fTrackBuffSize(19000),
282  fDoEventMixing(0),
284  fNzVtxBins (0),
285  fNCentBins (0),
286  fNRPBins (0),
287  fNOfPools(1),
288  fEventBuffer(0x0),
289  fEventInfo(0x0),
290  fProtonTracks(0x0),
291  fV0Tracks(0x0),
292  fProtonCutVarsArray(0x0),
293  fV0CutVarsArray(0x0)
294 {
295  //
296  // Constructor. Initialization of Inputs and Outputs
297  //
298  Info("AliAnalysisTaskSELc2pK0sfromAODtracks","Calling Constructor");
299 
300  DefineOutput(1,TList::Class()); //conters
301  DefineOutput(2,TList::Class());
302  DefineOutput(3,TList::Class()); //conters
303  DefineOutput(4,TTree::Class()); //My private output
304  DefineOutput(5,TTree::Class()); //My private output
305  DefineOutput(6,TTree::Class()); //My private output
306  DefineOutput(7,TTree::Class()); //My private output
307  DefineOutput(8,AliNormalizationCounter::Class());
308  DefineOutput(9,TTree::Class()); //My private output
309  DefineOutput(10,TTree::Class()); //My private output
310 }
311 
312 //___________________________________________________________________________
314  //
315  // destructor
316  //
317  Info("~AliAnalysisTaskSELc2pK0sfromAODtracks","Calling Destructor");
318 
319  if (fOutput) {
320  delete fOutput;
321  fOutput = 0;
322  }
323 
324  if (fOutputAll) {
325  delete fOutputAll;
326  fOutputAll = 0;
327  }
328 
329  if (fListCuts) {
330  delete fListCuts;
331  fListCuts = 0;
332  }
333 
334 
335  if (fAnalCuts) {
336  delete fAnalCuts;
337  fAnalCuts = 0;
338  }
339 
340  if (fVariablesTree) {
341  delete fVariablesTree;
342  fVariablesTree = 0;
343  }
344 
345  if (fProtonVariablesTree) {
346  delete fProtonVariablesTree;
348  }
349  if (fV0VariablesTree) {
350  delete fV0VariablesTree;
351  fV0VariablesTree = 0;
352  }
353  if (fMCVariablesTree) {
354  delete fMCVariablesTree;
355  fMCVariablesTree = 0;
356  }
358  delete fMCProtonVariablesTree;
360  }
361  if (fMCV0VariablesTree) {
362  delete fMCV0VariablesTree;
363  fMCV0VariablesTree = 0;
364  }
365 
366  if(fCounter){
367  delete fCounter;
368  fCounter = 0;
369  }
370 
371  if(fProtonTracks) fProtonTracks->Delete();
372  delete fProtonTracks;
373  if(fV0Tracks) fV0Tracks->Delete();
374  delete fV0Tracks;
376  delete fProtonCutVarsArray;
377  if(fV0CutVarsArray) fV0CutVarsArray->Delete();
378  delete fV0CutVarsArray;
379  if(fEventBuffer){
380  for(Int_t i=0; i<fNOfPools; i++) delete fEventBuffer[i];
381  delete fEventBuffer;
382  }
383  delete fEventInfo;
384 
385  if (fGTI)
386  delete[] fGTI;
387  fGTI=0;
388  if (fGTIndex)
389  delete[] fGTIndex;
390  fGTIndex=0;
391 
392 }
393 
394 //_________________________________________________
396  //
397  // Initialization
398  //
399  //
400 
401  fIsEventSelected=kFALSE;
402 
403  if (fDebug > 1) AliInfo("Init");
404 
405  fListCuts = new TList();
406  fListCuts->SetOwner();
407  fListCuts->SetName("ListCuts");
409  PostData(2,fListCuts);
410 
411  return;
412 }
413 
414 //_________________________________________________
416 {
417  //
418  // UserExec
419  //
420 
421  if (!fInputEvent) {
422  AliError("NO EVENT FOUND!");
423  return;
424  }
425  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
426  fCEvents->Fill(1);
428 
429  //------------------------------------------------
430  // First check if the event has proper vertex and B
431  //------------------------------------------------
432  fBzkG = (Double_t)aodEvent->GetMagneticField();
433  AliKFParticle::SetField(fBzkG);
434  if (TMath::Abs(fBzkG)<0.001) {
435  delete fV1;
436  return;
437  }
438  fCEvents->Fill(2);
439 
442 
443  //------------------------------------------------
444  // MC analysis setting
445  //------------------------------------------------
446  TClonesArray *mcArray = 0;
447  AliAODMCHeader *mcHeader=0;
448  if (fUseMCInfo) {
449  // MC array need for maching
450  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
451  if (!mcArray) {
452  AliError("Could not find Monte-Carlo in AOD");
453  return;
454  }
455  fCEvents->Fill(6); // in case of MC events
456 
457  // load MC header
458  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
459  if (!mcHeader) {
460  AliError("AliAnalysisTaskSELc2pK0sfromAODtracks::UserExec: MC header branch not found!\n");
461  return;
462  }
463  fCEvents->Fill(7); // in case of MC events
464 
465  Double_t zMCVertex = mcHeader->GetVtxZ();
466  if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()) {
467  AliDebug(2,Form("Event rejected: abs(zVtxMC)=%f > fAnalCuts->GetMaxVtxZ()=%f",zMCVertex,fAnalCuts->GetMaxVtxZ()));
468  return;
469  } else {
470  fCEvents->Fill(17); // in case of MC events
471  }
472  if ((TMath::Abs(zMCVertex) < fAnalCuts->GetMaxVtxZ()) && (!fAnalCuts->IsEventRejectedDuePhysicsSelection()) && (!fAnalCuts->IsEventRejectedDueToTrigger())) {
473  Bool_t selevt = MakeMCAnalysis(mcArray);
474  if(!selevt) return;
475  }
476  }
477 
478  //------------------------------------------------
479  // Event selection
480  //------------------------------------------------
481  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
482  if (!fVtx1) return;
483 
484  Double_t pos[3],cov[6];
485  fVtx1->GetXYZ(pos);
486  fVtx1->GetCovarianceMatrix(cov);
487  fV1 = new AliESDVertex(pos,cov,100.,100,fVtx1->GetName());
488  fVtxZ = pos[2];
489 
490  Bool_t fIsTriggerNotOK = fAnalCuts->IsEventRejectedDueToTrigger();
491  if(!fIsTriggerNotOK) fCEvents->Fill(3);
492  fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent); // better to initialize before CheckEventSelection call
493  if(!fIsEventSelected) {
494  delete fV1;
495  return;
496  }
497  fCEvents->Fill(4);
498 
499  fIsMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB)==(AliVEvent::kMB);
500  fIsSemi=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kSemiCentral)==(AliVEvent::kSemiCentral);
501  fIsCent=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kCentral)==(AliVEvent::kCentral);
502  fIsINT7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kINT7)==(AliVEvent::kINT7);
503  fIsEMC7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kEMC7)==(AliVEvent::kEMC7);
505  if(fIsMB) fHTrigger->Fill(1);
506  if(fIsSemi) fHTrigger->Fill(2);
507  if(fIsCent) fHTrigger->Fill(3);
508  if(fIsINT7) fHTrigger->Fill(4);
509  if(fIsEMC7) fHTrigger->Fill(5);
510  if(fIsMB|fIsSemi|fIsCent) fHTrigger->Fill(7);
511  if(fIsINT7|fIsEMC7) fHTrigger->Fill(8);
512  if(fIsMB&fIsSemi) fHTrigger->Fill(10);
513  if(fIsMB&fIsCent) fHTrigger->Fill(11);
514  if(fIsINT7&fIsEMC7) fHTrigger->Fill(12);
515 
516  if(fUseCentralityV0M){
517  AliCentrality *cent = aodEvent->GetCentrality();
518  fCentrality = cent->GetCentralityPercentile("V0M");
519  }else{
520  fCentrality = 1.;
521  }
522  if(fCentrality<0.||fCentrality>100.-0.0000001) {
523  delete fV1;
524  return;
525  }
526  fHCentrality->Fill(fCentrality);
527  fRunNumber = aodEvent->GetRunNumber();
528 
529  Int_t runnumber_offset = 0;
530  Int_t runnumber = aodEvent->GetRunNumber();
531  if(runnumber<=131000&&runnumber>=114000){
532  runnumber_offset = 114000;//lhc10bcde
533  }else if(runnumber<=196000&&runnumber>=195000){
534  runnumber_offset = 195000;//lhc13bc
535  }else if(runnumber<=170593&&runnumber>=167902){
536  runnumber_offset = 167902;//lhc11h
537  }
538  fHistonEvtvsRunNumber->Fill(runnumber-runnumber_offset,1.);
539 
540  //------------------------------------------------
541  // Check if the event has v0 candidate
542  //------------------------------------------------
543  //Int_t nv0 = aodEvent->GetNumberOfV0s();
544  fCEvents->Fill(5);
545 
546  //------------------------------------------------
547  // Main analysis done in this function
548  //------------------------------------------------
551  MakeAnalysis(aodEvent,mcArray);
552 
553 
554  PostData(1,fOutput);
555  PostData(3,fOutputAll);
556  PostData(4,fVariablesTree);
557  PostData(5,fProtonVariablesTree);
558  PostData(6,fV0VariablesTree);
559  PostData(7,fMCVariablesTree);
560  PostData(8,fCounter);
561  PostData(9,fMCProtonVariablesTree);
562  PostData(10,fMCV0VariablesTree);
563 
564  fIsEventSelected=kFALSE;
565 
566  delete fV1;
567  return;
568 }
569 
570 //________________________________________ terminate ___________________________
572 {
573  // The Terminate() function is the last function to be called during
574  // a query. It always runs on the client, it can be used to present
575  // the results graphically or save the results to file.
576 
577  //AliInfo("Terminate","");
578  AliAnalysisTaskSE::Terminate();
579 
580  fOutput = dynamic_cast<TList*> (GetOutputData(1));
581  if (!fOutput) {
582  AliError("fOutput not available");
583  return;
584  }
585 
586  fOutputAll = dynamic_cast<TList*> (GetOutputData(3));
587  if (!fOutputAll) {
588  AliError("fOutputAll not available");
589  return;
590  }
591 
592  return;
593 }
594 
595 //___________________________________________________________________________
597 {
598  //
599  // UserCreateOutputObject
600  //
601  //AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
602 
603  //------------------------------------------------
604  // output object setting
605  //------------------------------------------------
606  fOutput = new TList();
607  fOutput->SetOwner();
608  fOutput->SetName("chist0");
609  DefineGeneralHistograms(); // define general histograms
610  PostData(1,fOutput);
611 
612  fOutputAll = new TList();
613  fOutputAll->SetOwner();
614  fOutputAll->SetName("anahisto");
615  DefineAnalysisHistograms(); // define general histograms
616  PostData(3,fOutputAll);
617 
619  PostData(4,fVariablesTree);
620 
622  PostData(5,fProtonVariablesTree);
623 
625  PostData(6,fV0VariablesTree);
626 
628  PostData(7,fMCVariablesTree);
629 
631  PostData(9,fMCProtonVariablesTree);
632 
634  PostData(10,fMCV0VariablesTree);
635 
636  //Counter for Normalization
637  TString normName="NormalizationCounter";
638  AliAnalysisDataContainer *cont = GetOutputSlot(8)->GetContainer();
639  if(cont)normName=(TString)cont->GetName();
640  fCounter = new AliNormalizationCounter(normName.Data());
641  fCounter->Init();
642  PostData(8,fCounter);
643 
644  if(fDoEventMixing){
645  fProtonTracks = new TObjArray();
646  fProtonTracks->SetOwner();
647  fV0Tracks = new TObjArray();
648  fV0Tracks->SetOwner();
650  fProtonCutVarsArray->SetOwner();
651  fV0CutVarsArray = new TObjArray();
652  fV0CutVarsArray->SetOwner();
653 
655  fEventBuffer = new TTree*[fNOfPools];
656  for(Int_t i=0; i<fNOfPools; i++){
657  fEventBuffer[i]=new TTree(Form("EventBuffer_%d",i), "Temporary buffer for event mixing");
658  fEventBuffer[i]->Branch("zVertex", &fVtxZ);
659  fEventBuffer[i]->Branch("centrality", &fCentrality);
660  fEventBuffer[i]->Branch("reactionplane", &fReactionPlane);
661  fEventBuffer[i]->Branch("eventInfo", "TObjString",&fEventInfo);
662  fEventBuffer[i]->Branch("v1array", "TObjArray", &fV0Tracks);
663  fEventBuffer[i]->Branch("v1varsarray", "TObjArray", &fV0CutVarsArray);
664  }
665  }
666 
667  fGTI = new AliAODTrack *[fTrackBuffSize]; // Array of pointers
668  fGTIndex = new Int_t [fTrackBuffSize]; // Array of index
669 
670  return;
671 }
672 
673 //-------------------------------------------------------------------------------
675 (
676  AliAODEvent *aodEvent, TClonesArray *mcArray
677  )
678 {
679  //
680  // Main Analysis part
681  //
682  if(fDoEventMixing){
683  if(fProtonTracks) fProtonTracks->Delete();
684  if(fV0Tracks) fV0Tracks->Delete();
686  if(fV0CutVarsArray) fV0CutVarsArray->Delete();
687  }
688 
690  // ..and set it
691  for (Int_t iTrack=0;iTrack<aodEvent->GetNumberOfTracks();iTrack++){
692  // cast needed since the event now returns AliVTrack instead of AliAODTrack
693  AliAODTrack *track = dynamic_cast<AliAODTrack *>(aodEvent->GetTrack(iTrack));
694  if (!track) continue;
695 
696  // Store the reference of the global tracks
697  StoreGlobalTrackReference(track,iTrack);
698  }
699 
700 
701  Int_t nV0s= aodEvent->GetNumberOfV0s();
702  Int_t nTracks= aodEvent->GetNumberOfTracks();
703 
704  Int_t seleTrkFlags[nTracks];
705  Int_t nSeleTrks=0;
706  SelectTrack(aodEvent,nTracks,nSeleTrks,seleTrkFlags,mcArray);
707 
708  Bool_t seleV0Flags[nV0s];
709  Int_t nSeleV0=0;
710  SelectV0(aodEvent,nV0s,nSeleV0,seleV0Flags,mcArray);
711 
712  //------------------------------------------------
713  // V0 loop
714  //------------------------------------------------
715  for (Int_t iv0 = 0; iv0<nV0s; iv0++) {
716  if(!seleV0Flags[iv0]) continue;
717  AliAODv0 *v0 = aodEvent->GetV0(iv0);
718  if(!v0) continue;
719 
720  AliAODTrack *cptrack = (AliAODTrack*)(v0->GetDaughter(0));
721  AliAODTrack *cntrack = (AliAODTrack*)(v0->GetDaughter(1));
722 
723  //------------------------------------------------
724  // track loop
725  //------------------------------------------------
726  for (Int_t itrk = 0; itrk<nTracks; itrk++) {
727  if(seleTrkFlags[itrk]!=1) continue;
728  AliAODTrack *trk = (AliAODTrack*)aodEvent->GetTrack(itrk);
729  //if(trk->GetID()<0) continue;
730 
731  //TPC only track (BIT 7) does not have PID information
732  //In addition to that, TPC only tracks does not have good DCA resolution
733  //(according to femtoscopy code)
734  AliAODTrack *trkpid = 0;
735  if(fAnalCuts->GetProdAODFilterBit()==7){
736  trkpid = fGTI[-trk->GetID()-1];
737  }else{
738  trkpid = trk;
739  }
740 
741  Int_t cpid = cptrack->GetID();
742  Int_t cnid = cntrack->GetID();
743  Int_t lpid = trkpid->GetID();
744  if((cpid==lpid)||(cnid==lpid)) continue;
745 
746  if(!fAnalCuts->SelectWithRoughCuts(v0,trk)) continue;
747 
748  AliAODVertex *secVert = ReconstructSecondaryVertex(v0,trk,aodEvent);
749  if(!secVert) continue;
750 
751  AliAODRecoCascadeHF *lcobj = MakeCascadeHF(v0,trk,trkpid,aodEvent,secVert);
752  if(!lcobj) {
753  continue;
754  }
755 
756  FillROOTObjects(lcobj,v0,trk,trkpid,aodEvent,mcArray);
757 
758  lcobj->GetSecondaryVtx()->RemoveDaughters();
759  lcobj->UnsetOwnPrimaryVtx();
760  delete lcobj;lcobj=NULL;
761  delete secVert;
762  }
763  }
764 
765  if(fDoEventMixing){
766  fEventInfo->SetString(Form("Ev%d_esd%d_E%d_V%d",AliAnalysisManager::GetAnalysisManager()->GetNcalls(),((AliAODHeader*)aodEvent->GetHeader())->GetEventNumberESDFile(),fProtonTracks->GetEntries(),fV0Tracks->GetEntries()));
768  if(ind>=0 && ind<fNOfPools){
769  if(fEventBuffer[ind]->GetEntries() >= fNumberOfEventsForMixing){
771  if(fEventBuffer[ind]->GetEntries() >= 20*fNumberOfEventsForMixing){
772  ResetPool(ind);
773  }
774  }
775  fEventBuffer[ind]->Fill();
776  }
777  }
778 
779 }
780 
782 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillROOTObjects(AliAODRecoCascadeHF *lcobj, AliAODv0 *v0, AliAODTrack *trk, AliAODTrack *trkpid, AliAODEvent *aodEvent, TClonesArray *mcArray)
783 {
784  //
785  // Fill histograms or tree depending on fWriteVariableTree
786  //
787  if(!trk) return;
788  if(!trkpid) return;
789  if(!v0) return;
790 
791  Double_t mprPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
792  Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
793  Double_t mlamPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
794 
795  for(Int_t i=0;i<54;i++){
796  fCandidateVariables[i] = -9999.;
797  }
798 
799  Double_t pxp = trk->Px();
800  Double_t pyp = trk->Py();
801  Double_t pzp = trk->Pz();
802  Double_t momp = sqrt(pxp*pxp+pyp*pyp+pzp*pzp);
803  Double_t Ep = sqrt(momp*momp+mprPDG*mprPDG);
804 
805  Double_t pxv = v0->Px();
806  Double_t pyv = v0->Py();
807  Double_t pzv = v0->Pz();
808  Double_t momv = sqrt(pxv*pxv+pyv*pyv+pzv*pzv);
809  Double_t mv = v0->MassK0Short();
810  Double_t Ev = sqrt(momv*momv+mv*mv);
811 
812  Double_t cosoa = (pxp*pxv+pyp*pyv+pzp*pzv)/momp/momv;
813 
814  fCandidateVariables[ 0] = lcobj->InvMassLctoK0sP();
815  fCandidateVariables[ 1] = lcobj->Px();
816  fCandidateVariables[ 2] = lcobj->Py();
817  fCandidateVariables[ 3] = lcobj->Pz();
818  fCandidateVariables[ 4] = v0->MassK0Short();
819  fCandidateVariables[ 5] = lcobj->PxProng(0);
820  fCandidateVariables[ 6] = lcobj->PyProng(0);
821  fCandidateVariables[ 7] = lcobj->PzProng(0);
822  fCandidateVariables[ 8] = lcobj->PxProng(1);
823  fCandidateVariables[ 9] = lcobj->PyProng(1);
824  fCandidateVariables[10] = lcobj->PzProng(1);
825  fCandidateVariables[11] = fVtx1->GetX();
826  fCandidateVariables[12] = fVtx1->GetY();
827  fCandidateVariables[13] = fVtx1->GetZ();
829  fCandidateVariables[15] = lcobj->DecayLengthXY();
831 
832  Double_t nSigmaTPCpr=-9999.;
833  Double_t nSigmaTOFpr=-9999.;
834  Double_t probProton=-9999.;
835  if(fAnalCuts->GetIsUsePID())
836  {
837  nSigmaTPCpr = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTPC(trkpid,AliPID::kProton);
838  nSigmaTOFpr = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTOF(trkpid,AliPID::kProton);
840  probProton = fAnalCuts->GetProtonProbabilityTPCTOF(trk);
841  }
842  fCandidateVariables[17] = nSigmaTPCpr;
843  fCandidateVariables[18] = nSigmaTOFpr;
844  fCandidateVariables[19] = probProton;
845  }
846  fCandidateVariables[20] = 0;
847  fCandidateVariables[21] = lcobj->Getd0Prong(0);
848  fCandidateVariables[22] = lcobj->Getd0Prong(1);
849 
850  Double_t pospx = v0->MomPosX();
851  Double_t pospy = v0->MomPosY();
852  Double_t pospz = v0->MomPosZ();
853  Double_t pose = sqrt(pospx*pospx+pospy*pospy+pospz*pospz+0.000511*0.000511);
854  Double_t negpx = v0->MomNegX();
855  Double_t negpy = v0->MomNegY();
856  Double_t negpz = v0->MomNegZ();
857  Double_t nege = sqrt(negpx*negpx+negpy*negpy+negpz*negpz+0.000511*0.000511);
858  Double_t massPhoton = sqrt(pow(pose+nege,2)-pow(pospx+negpx,2)-pow(pospy+negpy,2)-pow(pospz+negpz,2));
859  fCandidateVariables[23] = v0->MassLambda();
860  fCandidateVariables[24] = v0->MassAntiLambda();
861  fCandidateVariables[25] = massPhoton;
862  fCandidateVariables[26] = v0->DcaPosToPrimVertex();
863  fCandidateVariables[27] = v0->DcaNegToPrimVertex();
864  fCandidateVariables[28] = v0->DcaV0ToPrimVertex();
865  fCandidateVariables[29] = v0->DcaV0Daughters();
866 
867  AliAODMCParticle *mclc = 0;
868  AliAODMCParticle *mcpr = 0;
869  AliAODMCParticle *mcv0 = 0;
870  Int_t mclablc = 0;
871  Int_t mcpdgpr_array[100];
872  Int_t mcpdgv0_array[100];
873  Int_t mclabelpr_array[100];
874  Int_t mclabelv0_array[100];
875  Int_t mcngen_pr=-9999;
876  Int_t mcngen_v0=-9999;
877  Bool_t islambdac = kFALSE;
878  Bool_t islambdac3body = kFALSE;
879  Double_t decayvertx_mc = -9999.;
880  Double_t decayverty_mc = -9999.;
881  Double_t decayvertz_mc = -9999.;
882  Double_t genvertx_mc = -9999.;
883  Double_t genverty_mc = -9999.;
884  Double_t genvertz_mc = -9999.;
885  if(fUseMCInfo && mcArray){
886  mclablc = MatchToMC(lcobj,mcArray,mcpdgpr_array, mcpdgv0_array,mclabelpr_array,mclabelv0_array,mcngen_pr,mcngen_v0);
887  if(mclablc>-1){
888  mclc = (AliAODMCParticle*) mcArray->At(mclablc);
889  if(mclabelpr_array[0]>=0)
890  mcpr = (AliAODMCParticle*) mcArray->At(mclabelpr_array[0]);
891  if(mclabelv0_array[0]>=0)
892  mcv0 = (AliAODMCParticle*) mcArray->At(mclabelv0_array[0]);
893  }
894  if(mclc){
895  Int_t pdgcode = mclc->GetPdgCode();
896  if(abs(pdgcode)==4122 && abs(mcpdgpr_array[1])==4122 && abs(mcpdgv0_array[1])==311 && abs(mcpdgv0_array[2])==4122 && mclc->GetNDaughters()==2){
897  islambdac = kTRUE;
898  }
899  if(abs(pdgcode)==4122 && abs(mcpdgpr_array[1])==2214 && abs(mcpdgpr_array[2])==4122 && abs(mcpdgv0_array[1])==311 && abs(mcpdgv0_array[2])==4122 && mclc->GetNDaughters()==2){
900  islambdac3body = kTRUE;
901  }
902  if(abs(pdgcode)==4122 && abs(mcpdgpr_array[1])==3222 && abs(mcpdgpr_array[2])==4122 && abs(mcpdgv0_array[1])==311 && abs(mcpdgv0_array[2])==4122 && mclc->GetNDaughters()==2){
903  islambdac3body = kTRUE;
904  }
905  if(abs(pdgcode)==4122 && abs(mcpdgpr_array[1])==3222 && abs(mcpdgpr_array[2])==3224 && abs(mcpdgpr_array[3])==4122 && abs(mcpdgv0_array[1])==311 && abs(mcpdgv0_array[2])==4122 && mclc->GetNDaughters()==2){
906  islambdac3body = kTRUE;
907  }
908  if(abs(pdgcode)==4122 && abs(mcpdgpr_array[1])==4122 && abs(mcpdgv0_array[1])==311 && abs(mcpdgv0_array[2])==313 && abs(mcpdgv0_array[3])==4122 && mclc->GetNDaughters()==2){
909  islambdac3body = kTRUE;
910  }
911  genvertx_mc = mclc->Xv();
912  genverty_mc = mclc->Yv();
913  genvertz_mc = mclc->Zv();
914  }
915  if(mcpr){
916  decayvertx_mc = mcpr->Xv();
917  decayverty_mc = mcpr->Yv();
918  decayvertz_mc = mcpr->Zv();
919  }
920  }
921  fCandidateVariables[30] = (Float_t) islambdac;
922 
923  Double_t LcPx = lcobj->Px();
924  Double_t LcPy = lcobj->Py();
925  Double_t LcPt = TMath::Sqrt(LcPx*LcPx+LcPy*LcPy);
926 
927  Double_t d0z0[2],covd0z0[3];
928  trk->PropagateToDCA(fVtx1,fBzkG,kVeryBig,d0z0,covd0z0);
929  Double_t x0 = fVtx1->GetX();
930  Double_t y0 = fVtx1->GetY();
931  Double_t px0 = LcPx/LcPt;
932  Double_t py0 = LcPy/LcPt;
933  Double_t tx[3];
934  trk->GetXYZ(tx);
935  Double_t x1 = tx[0];
936  Double_t y1 = tx[1];
937  trk->GetPxPyPz(tx);
938  Double_t px1 = tx[0];
939  Double_t py1 = tx[1];
940  Double_t pt1 = sqrt(px1*px1+py1*py1);
941  px1 = px1/pt1;
942  py1 = py1/pt1;
943 
944  Double_t dx = x0 - x1;
945  Double_t dy = y0 - y1;
946 
947  Double_t Delta = -px0*py1+py0*px1;
948  Double_t a0 = -9999.;
949  if(Delta!=0)
950  {
951  a0 = 1./Delta * (py1 * dx - px1 * dy);
952  }
953  Double_t neovertx = x0 + a0 * px0;
954  Double_t neoverty = y0 + a0 * py0;
955  Double_t z0 = fVtx1->GetZ();
956  Double_t neovertz = z0 + TMath::Abs(a0)*trk->Pz()/trk->Pt();
957 
958  fCandidateVariables[31] = neovertx;
959  fCandidateVariables[32] = neoverty;
960  fCandidateVariables[33] = neovertz;
961  fCandidateVariables[34] = (Float_t) islambdac3body;
962  fCandidateVariables[35] = decayvertx_mc;
963  fCandidateVariables[36] = decayverty_mc;
964  fCandidateVariables[37] = decayvertz_mc;
965  if(mclc){
966  fCandidateVariables[38] = mclc->Px();
967  fCandidateVariables[39] = mclc->Py();
968  fCandidateVariables[40] = mclc->Pz();
969  }
970  fCandidateVariables[41] = genvertx_mc;
971  fCandidateVariables[42] = genverty_mc;
972  fCandidateVariables[43] = genvertz_mc;
973 
974  Double_t x1_k0s = v0->DecayVertexV0X();
975  Double_t y1_k0s = v0->DecayVertexV0Y();
976  Double_t px1_k0s = v0->Px();
977  Double_t py1_k0s = v0->Py();
978  Double_t pt1_k0s = sqrt(px1_k0s*px1_k0s+py1_k0s*py1_k0s);
979  px1_k0s = px1_k0s/pt1_k0s;
980  py1_k0s = py1_k0s/pt1_k0s;
981 
982  Double_t dx_k0s = x0 - x1_k0s;
983  Double_t dy_k0s = y0 - y1_k0s;
984 
985  Double_t Delta_k0s = -px0*py1_k0s+py0*px1_k0s;
986  Double_t a0_k0s = -9999.;
987  if(Delta_k0s!=0)
988  {
989  a0_k0s = 1./Delta_k0s * (py1_k0s * dx_k0s - px1_k0s * dy_k0s);
990  }
991  Double_t neovertx_k0s = x0 + a0_k0s * px0;
992  Double_t neoverty_k0s = y0 + a0_k0s * py0;
993  fCandidateVariables[44] = neovertx_k0s;
994  fCandidateVariables[45] = neoverty_k0s;
995  fCandidateVariables[46] = lcobj->GetDCA();
996  if(mclc){
997  fCandidateVariables[47] = mclc->GetPdgCode();
998  fCandidateVariables[48] = mcpdgpr_array[0];
999  fCandidateVariables[49] = mcpdgpr_array[1];
1000  fCandidateVariables[50] = mcpdgpr_array[2];
1001  fCandidateVariables[51] = mcpdgv0_array[0];
1002  fCandidateVariables[52] = mcpdgv0_array[1];
1003  fCandidateVariables[53] = mcpdgv0_array[2];
1004  }
1005 
1006 
1007  if(fWriteVariableTree)
1008  fVariablesTree->Fill();
1009 
1010  Double_t cont_correlation[5];
1011  Double_t phi_trig = trk->Phi();
1012  Double_t phi_assoc = v0->Phi();
1013  Double_t eta_trig = trk->Eta();
1014  Double_t eta_assoc = v0->Eta();
1015  Double_t deltaphi = phi_assoc-phi_trig;
1016  if(deltaphi<-M_PI/2.) deltaphi += 2 * M_PI;
1017  if(deltaphi>3*M_PI/2.) deltaphi -= 2 * M_PI;
1018  cont_correlation[0] = deltaphi;
1019  cont_correlation[1] = eta_assoc-eta_trig;
1020  cont_correlation[2] = trk->Pt();
1021  cont_correlation[3] = v0->Pt();
1022  cont_correlation[4] = fCentrality;
1023  fHistoK0spCorrelation->Fill(cont_correlation);
1024 
1025  Double_t cont[3];
1026  cont[0] = lcobj->InvMassLctoK0sP();
1027  cont[1] = lcobj->Pt();
1028  cont[2] = fCentrality;
1029  if(cosoa>0.) fHistoLcK0SpMassCoarse->Fill(cont);
1030 
1032  {
1033  fHistoLcK0SpMass->Fill(cont);
1034 
1035  fHistod0Bach->Fill(lcobj->Getd0Prong(0));
1036  fHistod0V0->Fill(lcobj->Getd0Prong(1));
1037  fHistod0d0->Fill(lcobj->Getd0Prong(0)*lcobj->Getd0Prong(1));
1038  fHistoV0CosPA->Fill(lcobj->CosV0PointingAngle());
1039  fHistoProbProton->Fill(probProton);
1040  fHistoDecayLength->Fill(lcobj->DecayLengthXY()*(fAnalCuts->CalculateLcCosPAXY(lcobj)));
1041  fHistoK0SMass->Fill(v0->MassK0Short());
1042 
1043  if(fUseMCInfo){
1044  if(islambdac){
1045  fHistoLcK0SpMassMCS->Fill(cont);
1046  fHistoK0spCorrelationMCS->Fill(cont_correlation);
1047  fHistoResponseLcPt->Fill(mclc->Pt(),lcobj->Pt());
1048  }
1049  }
1050  }
1051  if(fUseMCInfo){
1052  if(islambdac3body){
1053  fHistoLcK0SpPi0MassMCS->Fill(cont);
1054  }
1055  }
1056 
1057  return;
1058 }
1060 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillMixROOTObjects(TLorentzVector *trkp, TLorentzVector *v0, TVector *prvars, TVector *v0vars)
1061 {
1065  if(!trkp) return;
1066  if(!v0) return;
1067 
1068  for(Int_t i=0;i<54;i++){
1069  fCandidateVariables[i] = -9999.;
1070  }
1071 
1072  Double_t pxp = trkp->Px();
1073  Double_t pyp = trkp->Py();
1074  Double_t pzp = trkp->Pz();
1075  Double_t momp = sqrt(pxp*pxp+pyp*pyp+pzp*pzp);
1076  Double_t Ep = sqrt(momp*momp+0.938272*0.938272);
1077 
1078  Double_t pxv = v0->Px();
1079  Double_t pyv = v0->Py();
1080  Double_t pzv = v0->Pz();
1081  Double_t momv = sqrt(pxv*pxv+pyv*pyv+pzv*pzv);
1082  Double_t mv = v0->M();
1083  Double_t Ev = sqrt(momv*momv+mv*mv);
1084 
1085  Double_t cosoa = (pxp*pxv+pyp*pyv+pzp*pzv)/momp/momv;
1086 
1087  Double_t pxsum = pxp + pxv;
1088  Double_t pysum = pyp + pyv;
1089  Double_t pzsum = pzp + pzv;
1090  Double_t Esum = Ep + Ev;
1091  Double_t mpk0s = sqrt(Esum*Esum-pxsum*pxsum-pysum*pysum-pzsum*pzsum);
1092 
1093 
1094  Double_t cont_correlation[5];
1095  Double_t phi_trig = trkp->Phi();
1096  Double_t phi_assoc = v0->Phi();
1097  Double_t eta_trig = trkp->Eta();
1098  Double_t eta_assoc = v0->Eta();
1099  Double_t deltaphi = phi_assoc-phi_trig;
1100  if(deltaphi<-M_PI/2.) deltaphi += 2 * M_PI;
1101  if(deltaphi>3*M_PI/2.) deltaphi -= 2 * M_PI;
1102  cont_correlation[0] = deltaphi;
1103  cont_correlation[1] = eta_assoc-eta_trig;
1104  cont_correlation[2] = trkp->Pt();
1105  cont_correlation[3] = v0->Pt();
1106  cont_correlation[4] = fCentrality;
1107  fHistoK0spCorrelationMix->Fill(cont_correlation);
1108 
1109  fCandidateVariables[ 0] = mpk0s;
1110  fCandidateVariables[ 1] = pxsum;
1111  fCandidateVariables[ 2] = pysum;
1112  fCandidateVariables[ 3] = pzsum;
1113  fCandidateVariables[ 4] = v0->M();
1114  fCandidateVariables[ 5] = trkp->Px();
1115  fCandidateVariables[ 6] = trkp->Py();
1116  fCandidateVariables[ 7] = trkp->Pz();
1117  fCandidateVariables[ 8] = v0->Px();
1118  fCandidateVariables[ 9] = v0->Py();
1119  fCandidateVariables[10] = v0->Pz();
1120  fCandidateVariables[20] = 1;
1121 
1122  Double_t LcPx = pxsum;
1123  Double_t LcPy = pysum;
1124  Double_t LcPt = TMath::Sqrt(LcPx*LcPx+LcPy*LcPy);
1125 
1126  Double_t x0 = (*prvars)[5];
1127  Double_t y0 = (*prvars)[6];
1128  Double_t px0 = LcPx/LcPt;
1129  Double_t py0 = LcPy/LcPt;
1130  Double_t tx[3];
1131  Double_t x1 = (*prvars)[1];
1132  Double_t y1 = (*prvars)[2];
1133  Double_t px1 = (*prvars)[3];
1134  Double_t py1 = (*prvars)[4];
1135 
1136  Double_t dx = x0 - x1;
1137  Double_t dy = y0 - y1;
1138 
1139  Double_t Delta = -px0*py1+py0*px1;
1140  Double_t a0 = -9999.;
1141  if(Delta!=0)
1142  {
1143  a0 = 1./Delta * (py1 * dx - px1 * dy);
1144  }
1145  Double_t neovertx = x0 + a0 * px0;
1146  Double_t neoverty = y0 + a0 * py0;
1147 
1148  fCandidateVariables[31] = neovertx;
1149  fCandidateVariables[32] = neoverty;
1150 
1151  if(fWriteVariableTree)
1152  fVariablesTree->Fill();
1153 
1154  Double_t rdhfcutvars[7];
1155  rdhfcutvars[0] = mpk0s;
1156  rdhfcutvars[1] = sqrt(pxsum*pxsum+pysum*pysum);
1157  rdhfcutvars[2] = (*prvars)[7];
1158  rdhfcutvars[3] = (*prvars)[8];
1159  rdhfcutvars[4] = a0;
1160  rdhfcutvars[5] = (*prvars)[0];
1161  rdhfcutvars[6] = (*v0vars)[0];
1162 
1163  if(fAnalCuts->IsSelected(trkp,v0,rdhfcutvars,AliRDHFCuts::kAll))
1164  {
1165  Double_t cont[3];
1166  cont[0] = mpk0s;
1167  cont[1] = sqrt(pxsum*pxsum+pysum*pysum);
1168  cont[2] = fCentrality;
1169  fHistoLcK0SpMassMixCoarse->Fill(cont);
1170  fHistoLcK0SpMassMix->Fill(cont);
1171  }
1172 
1173  return;
1174 }
1175 
1178 {
1179  //
1180  // Define tree variables
1181  //
1182 
1183  const char* nameoutput = GetOutputSlot(4)->GetContainer()->GetName();
1184  fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
1185  Int_t nVar = 54;
1187  TString * fCandidateVariableNames = new TString[nVar];
1188 
1189  fCandidateVariableNames[ 0]="InvMassLc2pK0s";
1190  fCandidateVariableNames[ 1]="LcPx";
1191  fCandidateVariableNames[ 2]="LcPy";
1192  fCandidateVariableNames[ 3]="LcPz";
1193  fCandidateVariableNames[ 4]="massK0Short";
1194  fCandidateVariableNames[ 5]="BachPx";
1195  fCandidateVariableNames[ 6]="BachPy";
1196  fCandidateVariableNames[ 7]="BachPz";
1197  fCandidateVariableNames[ 8]="V0Px";
1198  fCandidateVariableNames[ 9]="V0Py";
1199  fCandidateVariableNames[10]="V0Pz";
1200  fCandidateVariableNames[11]="PrimVertx";
1201  fCandidateVariableNames[12]="PrimVerty";
1202  fCandidateVariableNames[13]="PrimVertz";
1203  fCandidateVariableNames[14]="Centrality";
1204  fCandidateVariableNames[15]="DecayLengthXY";
1205  fCandidateVariableNames[16]="LcCosPAXY";
1206  fCandidateVariableNames[17]="nSigmaTPCpr";
1207  fCandidateVariableNames[18]="nSigmaTOFpr";
1208  fCandidateVariableNames[19]="probProton";
1209  fCandidateVariableNames[20]="Mixing";
1210  fCandidateVariableNames[21]="Bachd0";
1211  fCandidateVariableNames[22]="V0d0";
1212  fCandidateVariableNames[23]="massLambda";
1213  fCandidateVariableNames[24]="massAntiLambda";
1214  fCandidateVariableNames[25]="massPhoton";
1215  fCandidateVariableNames[26]="DcaPosToPrimVtx";
1216  fCandidateVariableNames[27]="DcaNegToPrimVtx";
1217  fCandidateVariableNames[28]="DcaV0ToPrimVtx";
1218  fCandidateVariableNames[29]="DcaV0Daughters";
1219  fCandidateVariableNames[30]="IsLambdac";
1220  fCandidateVariableNames[31]="SecVertX";
1221  fCandidateVariableNames[32]="SecVertY";
1222  fCandidateVariableNames[33]="SecVertZ";
1223  fCandidateVariableNames[34]="islambdac3body";
1224  fCandidateVariableNames[35]="SecVertXMC";
1225  fCandidateVariableNames[36]="SecVertYMC";
1226  fCandidateVariableNames[37]="SecVertZMC";
1227  fCandidateVariableNames[38]="LcPxMC";
1228  fCandidateVariableNames[39]="LcPyMC";
1229  fCandidateVariableNames[40]="LcPzMC";
1230  fCandidateVariableNames[41]="PrimVertXMC";
1231  fCandidateVariableNames[42]="PrimVertYMC";
1232  fCandidateVariableNames[43]="PrimVertZMC";
1233  fCandidateVariableNames[44]="SecVertXK0s";
1234  fCandidateVariableNames[45]="SecVertYK0s";
1235  fCandidateVariableNames[46]="DcaBachV0";
1236  fCandidateVariableNames[47]="MatchedPDG";
1237  fCandidateVariableNames[48]="ProtonPDG";
1238  fCandidateVariableNames[49]="MotherProtonPDG";
1239  fCandidateVariableNames[50]="GrMotherProtonPDG";
1240  fCandidateVariableNames[51]="V0PDG";
1241  fCandidateVariableNames[52]="MotherV0PDG";
1242  fCandidateVariableNames[53]="GrMotherV0PDG";
1243 
1244  for (Int_t ivar=0; ivar<nVar; ivar++) {
1245  fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
1246  }
1247 
1248  return;
1249 }
1250 
1253  //
1254  // This is to define general histograms
1255  //
1256 
1257  fCEvents = new TH1F("fCEvents","conter",18,-0.5,17.5);
1258  fCEvents->SetStats(kTRUE);
1259  fCEvents->GetXaxis()->SetBinLabel(1,"X1");
1260  fCEvents->GetXaxis()->SetBinLabel(2,"Analyzed events");
1261  fCEvents->GetXaxis()->SetBinLabel(3,"AliAODVertex exists");
1262  fCEvents->GetXaxis()->SetBinLabel(4,"TriggerOK");
1263  fCEvents->GetXaxis()->SetBinLabel(5,"IsEventSelected");
1264  fCEvents->GetXaxis()->SetBinLabel(6,"CascadesHF exists");
1265  fCEvents->GetXaxis()->SetBinLabel(7,"MCarray exists");
1266  fCEvents->GetXaxis()->SetBinLabel(8,"MCheader exists");
1267  fCEvents->GetXaxis()->SetBinLabel(9,"triggerClass!=CINT1");
1268  fCEvents->GetXaxis()->SetBinLabel(10,"triggerMask!=kAnyINT");
1269  fCEvents->GetXaxis()->SetBinLabel(11,"triggerMask!=kAny");
1270  fCEvents->GetXaxis()->SetBinLabel(12,"vtxTitle.Contains(Z)");
1271  fCEvents->GetXaxis()->SetBinLabel(13,"vtxTitle.Contains(3D)");
1272  fCEvents->GetXaxis()->SetBinLabel(14,"vtxTitle.Doesn'tContain(Z-3D)");
1273  fCEvents->GetXaxis()->SetBinLabel(15,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
1274  fCEvents->GetXaxis()->SetBinLabel(16,"!IsEventSelected");
1275  fCEvents->GetXaxis()->SetBinLabel(17,"triggerMask!=kAnyINT || triggerClass!=CINT1");
1276  fCEvents->GetXaxis()->SetBinLabel(18,Form("zVtxMC<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
1277  //fCEvents->GetXaxis()->SetTitle("");
1278  fCEvents->GetYaxis()->SetTitle("counts");
1279 
1280  fHTrigger = new TH1F("fHTrigger","counter",18,-0.5,17.5);
1281  fHTrigger->SetStats(kTRUE);
1282  fHTrigger->GetXaxis()->SetBinLabel(1,"X1");
1283  fHTrigger->GetXaxis()->SetBinLabel(2,"kMB");
1284  fHTrigger->GetXaxis()->SetBinLabel(3,"kSemiCentral");
1285  fHTrigger->GetXaxis()->SetBinLabel(4,"kCentral");
1286  fHTrigger->GetXaxis()->SetBinLabel(5,"kINT7");
1287  fHTrigger->GetXaxis()->SetBinLabel(6,"kEMC7");
1288  //fHTrigger->GetXaxis()->SetBinLabel(7,"Space");
1289  fHTrigger->GetXaxis()->SetBinLabel(8,"kMB|kSemiCentral|kCentral");
1290  fHTrigger->GetXaxis()->SetBinLabel(9,"kINT7|kEMC7");
1291  fHTrigger->GetXaxis()->SetBinLabel(11,"kMB&kSemiCentral");
1292  fHTrigger->GetXaxis()->SetBinLabel(12,"kMB&kCentral");
1293  fHTrigger->GetXaxis()->SetBinLabel(13,"kINT7&kEMC7");
1294 
1295  fHCentrality = new TH1F("fHCentrality","conter",100,0.,100.);
1296  fHReactionPlane = new TH1F("fHReactionPlane","conter",100,-1.6,3.2);
1297 
1298  fOutput->Add(fCEvents);
1299  fOutput->Add(fHTrigger);
1300  fOutput->Add(fHCentrality);
1301  fOutput->Add(fHReactionPlane);
1302 
1303  return;
1304 }
1305 //__________________________________________________________________________
1307 {
1308  //
1309  // Define analyis histograms
1310  //
1311 
1312  //------------------------------------------------
1313  // Basic histogram
1314  //------------------------------------------------
1315  Int_t bins_base[3]= {200 ,20 ,10};
1316  Double_t xmin_base[3]={2.286-0.5,0 ,0.00};
1317  Double_t xmax_base[3]={2.286+0.5,20. ,100};
1318  fHistoLcK0SpMass = new THnSparseF("fHistoLcK0SpMass","",3,bins_base,xmin_base,xmax_base);
1320  fHistoLcK0SpMassMix = new THnSparseF("fHistoLcK0SpMassMix","",3,bins_base,xmin_base,xmax_base);
1322  fHistoLcK0SpMassMCS = new THnSparseF("fHistoLcK0SpMassMCS","",3,bins_base,xmin_base,xmax_base);
1324  fHistoLcK0SpPi0MassMCS = new THnSparseF("fHistoLcK0SpPi0MassMCS","",3,bins_base,xmin_base,xmax_base);
1326  fHistoLcKPluspMass = new THnSparseF("fHistoLcKPluspMass","",3,bins_base,xmin_base,xmax_base);
1328  fHistoLcKMinuspMass = new THnSparseF("fHistoLcKMinuspMass","",3,bins_base,xmin_base,xmax_base);
1330  fHistoLcKPluspMassMix = new THnSparseF("fHistoLcKPluspMassMix","",3,bins_base,xmin_base,xmax_base);
1332  fHistoLcKMinuspMassMix = new THnSparseF("fHistoLcKMinuspMassMix","",3,bins_base,xmin_base,xmax_base);
1334 
1335  Int_t bins_lcmcgen[3]= {100 ,20 ,10};
1336  Double_t xmin_lcmcgen[3]={0.,-1.0 ,0.0};
1337  Double_t xmax_lcmcgen[3]={20.,1.0 ,100};
1338  fHistoLcMCGen = new THnSparseF("fHistoLcMCGen","",3,bins_lcmcgen,xmin_lcmcgen,xmax_lcmcgen);
1339  fOutputAll->Add(fHistoLcMCGen);
1340 
1341  Int_t bins_coarse[3]= {160 ,20 ,10};
1342  Double_t xmin_coarse[3]={1.,0 ,0.00};
1343  Double_t xmax_coarse[3]={5.8,20. ,100};
1344  fHistoLcK0SpMassCoarse = new THnSparseF("fHistoLcK0SpMassCoarse","",3,bins_coarse,xmin_coarse,xmax_coarse);
1346  fHistoLcK0SpMassMixCoarse = new THnSparseF("fHistoLcK0SpMassMixCoarse","",3,bins_coarse,xmin_coarse,xmax_coarse);
1348 
1349  Int_t bins_correlation[5]= {40 ,20,20,20,10};
1350  Double_t xmin_correlation[5]={-M_PI/2.,-2.4,0.,0.,0.0};
1351  Double_t xmax_correlation[5]={1.5*M_PI,2.4,10.,10.,100};
1352  fHistoK0spCorrelation = new THnSparseF("fHistoK0spCorrelation","",5,bins_correlation,xmin_correlation,xmax_correlation);
1354  fHistoK0spCorrelationMix = new THnSparseF("fHistoK0spCorrelationMix","",5,bins_correlation,xmin_correlation,xmax_correlation);
1356  fHistoK0spCorrelationMCS = new THnSparseF("fHistoK0spCorrelationMCS","",5,bins_correlation,xmin_correlation,xmax_correlation);
1358 
1359 
1360  //------------------------------------------------
1361  // checking histograms
1362  //------------------------------------------------
1363  fHistoBachPt = new TH2F("fHistoBachPt","Bachelor p_{T}",100,0.0,10.0,20,-1.,1.);
1364  fOutputAll->Add(fHistoBachPt);
1365  fHistoBachPtMCS = new TH2F("fHistoBachPtMCS","Bachelor p_{T}",100,0.0,10.0,20,-1.,1.);
1367  fHistoBachPtMCGen = new TH2F("fHistoBachPtMCGen","Bachelor p_{T}",100,0.0,10.0,20,-1.,1.);
1369 
1370  fHistoKaonPt = new TH2F("fHistoKaonPt","Kaon p_{T}",100,0.0,10.0,20,-1.,1.);
1371  fOutputAll->Add(fHistoKaonPt);
1372  fHistoKaonPtMCS = new TH2F("fHistoKaonPtMCS","Kaon p_{T}",100,0.0,10.0,20,-1.,1.);
1374  fHistoKaonPtMCGen = new TH2F("fHistoKaonPtMCGen","Kaon p_{T}",100,0.0,10.0,20,-1.,1.);
1376 
1377  fHistoK0sMassvsPt=new TH3F("fHistoK0sMassvsPt","K0s mass",100,0.497-0.05,0.497+0.05,20,0.,10.,20,-1.,1.);
1379  fHistoK0sMassvsPtMCS=new TH3F("fHistoK0sMassvsPtMCS","K0s mass",100,0.497-0.05,0.497+0.05,20,0.,10.,20,-1.,1.);
1381  fHistoK0sMassvsPtMCGen=new TH3F("fHistoK0sMassvsPtMCGen","K0s mass",100,0.497-0.05,0.497+0.05,20,0.,10.,20,-1.,1.);
1383 
1384  fHistod0Bach = new TH1F("fHistod0Bach","Bachelor d_{0}",100,-0.5,0.5);
1385  fOutputAll->Add(fHistod0Bach);
1386  fHistod0V0 = new TH1F("fHistod0V0","V_{0} d_{0}",100,-0.5,0.5);
1387  fOutputAll->Add(fHistod0V0);
1388  fHistod0d0 = new TH1F("fHistod0d0","d_{0} d_{0}",100,-0.5,0.5);
1389  fOutputAll->Add(fHistod0d0);
1390  fHistoV0CosPA=new TH1F("fHistoV0CosPA","V0->Second vertex cospa",100,-1.,1.0);
1391  fOutputAll->Add(fHistoV0CosPA);
1392  fHistoProbProton=new TH1F("fHistoProbProton","ProbProton",100,0.,1.0);
1394  fHistoDecayLength=new TH1F("fHistoDecayLength","Decay Length",100,-0.1,0.1);
1396  fHistoK0SMass=new TH1F("fHistoK0SMass","K0S mass",100,0.497-0.05,0.497+0.05);
1397  fOutputAll->Add(fHistoK0SMass);
1398 
1399  fHistonEvtvsRunNumber=new TH1F("fHistonEvtvsRunNumber","",20000,-0.5,19999.5);
1401  fHistonProtonvsRunNumber=new TH1F("fHistonProtonvsRunNumber","",20000,-0.5,19999.5);
1403  fHistonK0svsRunNumber=new TH1F("fHistonK0svsRunNumber","",20000,-0.5,19999.5);
1405 
1406  fHistoMassTagV0Min=new TH1F("fHistoMassTagV0Min","",1500,0,1.5);
1408  fHistoMassTagV0SameSignMin=new TH1F("fHistoMassTagV0SameSignMin","",1500,0,1.5);
1410 
1411  fHistoResponseLcPt = new TH2D("fHistoResponseLcPt","",100,0.,20.,100,0.,20.);
1413  fHistoResponseLcPt1 = new TH2D("fHistoResponseLcPt1","",100,0.,20.,100,0.,20.);
1415  fHistoResponseLcPt2 = new TH2D("fHistoResponseLcPt2","",100,0.,20.,100,0.,20.);
1417 
1418  return;
1419 }
1420 
1421 //________________________________________________________________________
1422 AliAODRecoCascadeHF* AliAnalysisTaskSELc2pK0sfromAODtracks::MakeCascadeHF(AliAODv0 *v0, AliAODTrack *part, AliAODTrack *partpid, AliAODEvent * aod, AliAODVertex *secVert)
1423 {
1424  //
1425  // Create AliAODRecoCascadeHF object from the argument
1426  //
1427 
1428  if(!v0) return 0x0;
1429  if(!part) return 0x0;
1430  if(!partpid) return 0x0;
1431  if(!aod) return 0x0;
1432 
1433  //------------------------------------------------
1434  // PrimaryVertex
1435  //------------------------------------------------
1436  AliAODVertex *primVertexAOD;
1437  Bool_t unsetvtx = kFALSE;
1439  primVertexAOD = CallPrimaryVertex(v0,part,aod);
1440  if(!primVertexAOD){
1441  primVertexAOD = fVtx1;
1442  }else{
1443  unsetvtx = kTRUE;
1444  }
1445  }else{
1446  primVertexAOD = fVtx1;
1447  }
1448  if(!primVertexAOD) return 0x0;
1449  Double_t posprim[3]; primVertexAOD->GetXYZ(posprim);
1450 
1451  //------------------------------------------------
1452  // DCA between tracks
1453  //------------------------------------------------
1454  AliESDtrack *esdtrack = new AliESDtrack((AliVTrack*)partpid);
1455 
1456  AliNeutralTrackParam *trackV0=NULL;
1457  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
1458  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
1459 
1460  Double_t xdummy, ydummy;
1461  Double_t dca = esdtrack->GetDCA(trackV0,fBzkG,xdummy,ydummy);
1462 
1463 
1464  //------------------------------------------------
1465  // Propagate all tracks to the secondary vertex and calculate momentum there
1466  //------------------------------------------------
1467 
1468  Double_t d0z0bach[2],covd0z0bach[3];
1469 // if(sqrt(pow(secVert->GetX(),2)+pow(secVert->GetY(),2))<1.){
1470 // part->PropagateToDCA(secVert,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
1471 // trackV0->PropagateToDCA(secVert,fBzkG,kVeryBig);
1472 // }else{
1473 // part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
1474 // trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig);
1475 // }
1476  part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
1477  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig);
1478  Double_t momv0_new[3]={-9999,-9999,-9999.};
1479  trackV0->GetPxPyPz(momv0_new);
1480 
1481  Double_t px[2],py[2],pz[2];
1482  px[0] = part->Px(); py[0] = part->Py(); pz[0] = part->Pz();
1483  px[1] = momv0_new[0]; py[1] = momv0_new[1]; pz[1] = momv0_new[2];
1484 
1485  //------------------------------------------------
1486  // d0
1487  //------------------------------------------------
1488  Double_t d0[3],d0err[3];
1489 
1490  part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
1491  d0[0]= d0z0bach[0];
1492  d0err[0] = TMath::Sqrt(covd0z0bach[0]);
1493 
1494  Double_t d0z0v0[2],covd0z0v0[3];
1495  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0v0,covd0z0v0);
1496  d0[1]= d0z0v0[0];
1497  d0err[1] = TMath::Sqrt(covd0z0v0[0]);
1498 
1499  //------------------------------------------------
1500  // Create AliAODRecoCascadeHF
1501  //------------------------------------------------
1502  Short_t charge = part->Charge();
1503  AliAODRecoCascadeHF *theCascade = new AliAODRecoCascadeHF(secVert,charge,px,py,pz,d0,d0err,dca);
1504  if(!theCascade)
1505  {
1506  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1507  if(esdtrack) delete esdtrack;
1508  if(trackV0) delete trackV0;
1509  return 0x0;
1510  }
1511  theCascade->SetOwnPrimaryVtx(primVertexAOD);
1512  UShort_t id[2]={(UShort_t)part->GetID(),(UShort_t)trackV0->GetID()};
1513  theCascade->SetProngIDs(2,id);
1514 
1515  theCascade->GetSecondaryVtx()->AddDaughter(partpid);
1516  theCascade->GetSecondaryVtx()->AddDaughter(v0);
1517 
1518  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1519  if(esdtrack) delete esdtrack;
1520  if(trackV0) delete trackV0;
1521 
1522  return theCascade;
1523 }
1524 
1525 //________________________________________________________________________
1526 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::CallPrimaryVertex(AliAODv0 *v0, AliAODTrack *trk, AliAODEvent* aod)
1527 {
1528  //
1529  // Make an array of tracks which should not be used in primary vertex calculation and
1530  // Call PrimaryVertex function
1531  //
1532 
1533  TObjArray *TrackArray = new TObjArray(3);
1534 
1535  AliESDtrack *cptrk1 = new AliESDtrack((AliVTrack*)trk);
1536  TrackArray->AddAt(cptrk1,0);
1537 
1538  AliESDtrack *cascptrack = new AliESDtrack((AliVTrack*)v0->GetDaughter(0));
1539  TrackArray->AddAt(cascptrack,1);
1540  AliESDtrack *cascntrack = new AliESDtrack((AliVTrack*)v0->GetDaughter(1));
1541  TrackArray->AddAt(cascntrack,2);
1542 
1543  AliAODVertex *newvert = PrimaryVertex(TrackArray,aod);
1544 
1545  for(Int_t i=0;i<3;i++)
1546  {
1547  AliESDtrack *tesd = (AliESDtrack*)TrackArray->UncheckedAt(i);
1548  delete tesd;
1549  }
1550  TrackArray->Clear();
1551  delete TrackArray;
1552 
1553  return newvert;
1554 }
1555 
1556 //________________________________________________________________________
1558  AliVEvent *event)
1559 {
1560  //
1561  //Used only for pp
1562  //copied from AliAnalysisVertexingHF (except for the following 3 lines)
1563  //
1564 
1565  Bool_t fRecoPrimVtxSkippingTrks = kTRUE;
1566  Bool_t fRmTrksFromPrimVtx = kFALSE;
1567 
1568  AliESDVertex *vertexESD = 0;
1569  AliAODVertex *vertexAOD = 0;
1570 
1571  //vertexESD = new AliESDVertex(*fV1);
1572 
1573 
1574  if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
1575  // primary vertex from the input event
1576 
1577  vertexESD = new AliESDVertex(*fV1);
1578 
1579  } else {
1580  // primary vertex specific to this candidate
1581 
1582  Int_t nTrks = trkArray->GetEntriesFast();
1583  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1584 
1585  if(fRecoPrimVtxSkippingTrks) {
1586  // recalculating the vertex
1587 
1588  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
1589  Float_t diamondcovxy[3];
1590  event->GetDiamondCovXY(diamondcovxy);
1591  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1592  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1593  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1594  vertexer->SetVtxStart(diamond);
1595  delete diamond; diamond=NULL;
1596  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
1597  vertexer->SetOnlyFitter();
1598  }
1599  Int_t skipped[1000];
1600  Int_t nTrksToSkip=0,id;
1601  AliExternalTrackParam *t = 0;
1602  for(Int_t i=0; i<nTrks; i++) {
1603  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
1604  id = (Int_t)t->GetID();
1605  if(id<0) continue;
1606  skipped[nTrksToSkip++] = id;
1607  }
1608  // TEMPORARY FIX
1609  // For AOD, skip also tracks without covariance matrix
1610  Double_t covtest[21];
1611  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
1612  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
1613  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
1614  id = (Int_t)vtrack->GetID();
1615  if(id<0) continue;
1616  skipped[nTrksToSkip++] = id;
1617  }
1618  }
1619  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
1620  //
1621  vertexer->SetSkipTracks(nTrksToSkip,skipped);
1622  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1623 
1624  } else if(fRmTrksFromPrimVtx && nTrks>0) {
1625  // removing the prongs tracks
1626 
1627  TObjArray rmArray(nTrks);
1628  UShort_t *rmId = new UShort_t[nTrks];
1629  AliESDtrack *esdTrack = 0;
1630  AliESDtrack *t = 0;
1631  for(Int_t i=0; i<nTrks; i++) {
1632  t = (AliESDtrack*)trkArray->UncheckedAt(i);
1633  esdTrack = new AliESDtrack(*t);
1634  rmArray.AddLast(esdTrack);
1635  if(esdTrack->GetID()>=0) {
1636  rmId[i]=(UShort_t)esdTrack->GetID();
1637  } else {
1638  rmId[i]=9999;
1639  }
1640  }
1641  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
1642  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
1643  delete [] rmId; rmId=NULL;
1644  rmArray.Delete();
1645 
1646  }
1647 
1648  delete vertexer; vertexer=NULL;
1649  if(!vertexESD) return vertexAOD;
1650  if(vertexESD->GetNContributors()<=0) {
1651  //AliDebug(2,"vertexing failed");
1652  delete vertexESD; vertexESD=NULL;
1653  return vertexAOD;
1654  }
1655 
1656 
1657  }
1658 
1659  // convert to AliAODVertex
1660  Double_t pos[3],cov[6],chi2perNDF;
1661  vertexESD->GetXYZ(pos); // position
1662  vertexESD->GetCovMatrix(cov); //covariance matrix
1663  chi2perNDF = vertexESD->GetChi2toNDF();
1664  delete vertexESD; vertexESD=NULL;
1665 
1666  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1667 
1668  return vertexAOD;
1669 }
1670 
1671 //________________________________________________________________________
1673 {
1674  //
1675  // Reconstruct secondary vertex from trkArray (Copied from AliAnalysisVertexingHF)
1676  //
1677  //------------------------------------------------
1678  // PrimaryVertex
1679  //------------------------------------------------
1680  AliAODVertex *primVertexAOD;
1681  Bool_t unsetvtx = kFALSE;
1683  primVertexAOD = CallPrimaryVertex(v0,part,aod);
1684  if(!primVertexAOD){
1685  primVertexAOD = fVtx1;
1686  }else{
1687  unsetvtx = kTRUE;
1688  }
1689  }else{
1690  primVertexAOD = fVtx1;
1691  }
1692  if(!primVertexAOD) return 0x0;
1693 
1694  //------------------------------------------------
1695  // Secondary vertex
1696  //------------------------------------------------
1697 
1698  Double_t LcPx = part->Px()+v0->Px();
1699  Double_t LcPy = part->Py()+v0->Py();
1700  Double_t LcPt = TMath::Sqrt(LcPx*LcPx+LcPy*LcPy);
1701 
1702  Double_t d0z0[2],covd0z0[3];
1703  part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1704  Double_t x0 = primVertexAOD->GetX();
1705  Double_t y0 = primVertexAOD->GetY();
1706  Double_t px0 = LcPx/LcPt;
1707  Double_t py0 = LcPy/LcPt;
1708  Double_t tx[3];
1709  part->GetXYZ(tx);
1710  Double_t x1 = tx[0];
1711  Double_t y1 = tx[1];
1712  part->GetPxPyPz(tx);
1713  Double_t px1 = tx[0];
1714  Double_t py1 = tx[1];
1715  Double_t pt1 = sqrt(px1*px1+py1*py1);
1716  px1 = px1/pt1;
1717  py1 = py1/pt1;
1718 
1719  Double_t dx = x0 - x1;
1720  Double_t dy = y0 - y1;
1721 
1722  Double_t Delta = -px0*py1+py0*px1;
1723  Double_t a0 = -9999.;
1724  if(Delta!=0)
1725  {
1726  a0 = 1./Delta * (py1 * dx - px1 * dy);
1727  }
1728  Double_t neovertx = x0 + a0 * px0;
1729  Double_t neoverty = y0 + a0 * py0;
1730  Double_t z0 = primVertexAOD->GetZ();
1731  Double_t neovertz = z0 + TMath::Abs(a0)*part->Pz()/part->Pt();
1732 
1733  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1734 
1735  Double_t pos[3],cov[6],chi2perNDF;
1736  pos[0]=neovertx;
1737  pos[1]=neoverty;
1738  pos[2]=neovertz;
1739  cov[0]=0.0;
1740  cov[1]=0.0;
1741  cov[2]=0.0;
1742  cov[3]=0.0;
1743  cov[4]=0.0;
1744  cov[5]=0.0;
1745  chi2perNDF=0.0;
1746  AliAODVertex *secVert = new AliAODVertex(pos,cov,chi2perNDF);
1747  if(!secVert){
1748  return 0x0;
1749  }
1750  return secVert;
1751 }
1752 //________________________________________________________________________
1753 //AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::ReconstructSecondaryVertex(AliAODv0 *v0, AliAODTrack *part, AliAODEvent * aod)
1754 //{
1755 // //
1756 // // Reconstruct secondary vertex from trkArray (Copied from AliAnalysisVertexingHF)
1757 // // Currently only returns Primary vertex (can we reconstruct secondary vertex from p - v0)
1758 // //
1759 //
1760 // AliAODVertex *primVertexAOD;
1761 // Bool_t unsetvtx = kFALSE;
1762 // if(fReconstructPrimVert){
1763 // primVertexAOD = CallPrimaryVertex(v0,part,aod);
1764 // if(!primVertexAOD){
1765 // primVertexAOD = fVtx1;
1766 // }else{
1767 // unsetvtx = kTRUE;
1768 // }
1769 // }else{
1770 // primVertexAOD = fVtx1;
1771 // }
1772 // if(!primVertexAOD) return 0x0;
1773 //
1774 // AliESDVertex * vertexESD = new AliESDVertex(*fV1);
1775 //
1776 // Double_t pos[3],cov[6],chi2perNDF;
1777 // vertexESD->GetXYZ(pos); // position
1778 // vertexESD->GetCovMatrix(cov); //covariance matrix
1779 // chi2perNDF = vertexESD->GetChi2toNDF();
1780 // delete vertexESD; vertexESD=NULL;
1781 //
1782 // AliAODVertex *secVert = new AliAODVertex(pos,cov,chi2perNDF);
1783 //
1784 // return secVert;
1785 //}
1786 //________________________________________________________________________
1787 void AliAnalysisTaskSELc2pK0sfromAODtracks::SelectTrack( const AliVEvent *event, Int_t trkEntries, Int_t &nSeleTrks,Int_t *seleFlags, TClonesArray *mcArray)
1788 {
1789  //
1790  // Select good tracks using fAnalCuts (AliRDHFCuts object) and return the array of their ids
1791  //
1792 
1793  if(trkEntries==0) return;
1794 
1795  nSeleTrks=0;
1796  for(Int_t i=0; i<trkEntries; i++) {
1797  seleFlags[i] =0;
1798 
1799  AliVTrack *track;
1800  track = (AliVTrack*)event->GetTrack(i);
1801 
1802  //if(track->GetID()<0) continue;
1803  Double_t covtest[21];
1804  if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
1805  if(!fAnalCuts) continue;
1806 
1807  AliAODTrack *aodt = (AliAODTrack*)track;
1808 
1810  Int_t filterbit = fAnalCuts->GetProdAODFilterBit();
1811  if(filterbit==7){
1812  if(!aodt->TestFilterBit(BIT(filterbit))) continue;
1813  }else{
1814  if(!aodt->TestFilterMask(BIT(filterbit))) continue;
1815  }
1816  }
1817 
1818  AliAODTrack *aodtpid = 0;
1819  if(fAnalCuts->GetProdAODFilterBit()==7){
1820  aodtpid = fGTI[-aodt->GetID()-1];
1821  }else{
1822  aodtpid = aodt;
1823  }
1824 
1825  if(fAnalCuts->SingleTrkCuts(aodt,aodtpid,fVtx1)){
1826  seleFlags[i]=1;
1827  nSeleTrks++;
1828  FillProtonROOTObjects(aodt,mcArray);
1829 
1830  Double_t minmass;
1831  Bool_t isv0 = fAnalCuts->TagV0(aodt,(AliAODEvent*)event,trkEntries,minmass);
1832  fHistoMassTagV0Min->Fill(minmass);
1833  if(isv0) seleFlags[i] = 0;
1834 
1835 // Double_t minmasslike;
1836 // fAnalCuts->TagV0SameSign(aodt,(AliAODEvent*)event,trkEntries,minmasslike);
1837 // fHistoMassTagV0SameSignMin->Fill(minmasslike);
1838  }
1839  } // end loop on tracks
1840 }
1841 //________________________________________________________________________
1842 void AliAnalysisTaskSELc2pK0sfromAODtracks::SelectV0( const AliVEvent *event,Int_t nV0s,Int_t &nSeleV0, Bool_t *seleV0Flags, TClonesArray *mcArray)
1843 {
1844  //
1845  // Select good V0 using fAnalCuts (AliRDHFCuts object) and return the array of their ids
1846  //
1847 
1848  nSeleV0 = 0;
1849  for(Int_t iv0=0;iv0<nV0s;iv0++)
1850  {
1851  seleV0Flags[iv0] = kFALSE;
1852  AliAODv0 *v0 = ((AliAODEvent*)event)->GetV0(iv0);
1853 
1854  if(!fAnalCuts) continue;
1855  if(fAnalCuts->SingleV0Cuts(v0,fVtx1)){
1856  seleV0Flags[iv0] = kTRUE;
1857  nSeleV0++;
1858 
1859  FillV0ROOTObjects(v0, mcArray);
1860  }
1861  }
1862 }
1863 
1866 {
1867  //
1869  //
1870 
1871  const char* nameoutput = GetOutputSlot(5)->GetContainer()->GetName();
1872  fProtonVariablesTree = new TTree(nameoutput,"proton variables tree");
1873  Int_t nVar = 3;
1875  TString * fCandidateVariableNames = new TString[nVar];
1876 
1877  fCandidateVariableNames[ 0]="PrPx";
1878  fCandidateVariableNames[ 1]="PrPy";
1879  fCandidateVariableNames[ 2]="PrPz";
1880 
1881  for (Int_t ivar=0; ivar<nVar; ivar++) {
1882  fProtonVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateProtonVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
1883  }
1884 
1885  return;
1886 }
1887 //________________________________________________________________________
1888 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillProtonROOTObjects(AliAODTrack *trk, TClonesArray *mcArray)
1889 {
1890  //
1892  //
1893 
1894  if(!trk) return;
1895 
1896  TLorentzVector vpr;
1897  vpr.SetXYZM(trk->Px(),trk->Py(),trk->Pz(),0.938272);
1898  fHistoBachPt->Fill(vpr.Pt(),vpr.Rapidity());
1899  if(fDoEventMixing){
1900  fProtonTracks->AddLast(new TLorentzVector(trk->Px(),trk->Py(),trk->Pz(),trk->Charge()));
1901 
1902  Double_t d0z0[2],covd0z0[3];
1903  trk->PropagateToDCA(fVtx1,fBzkG,kVeryBig,d0z0,covd0z0);
1904  Double_t tx[3];
1905  trk->GetXYZ(tx);
1906  Double_t x1 = tx[0];
1907  Double_t y1 = tx[1];
1908  trk->GetPxPyPz(tx);
1909  Double_t px1 = tx[0];
1910  Double_t py1 = tx[1];
1911  Double_t pt1 = sqrt(px1*px1+py1*py1);
1912  px1 = px1/pt1;
1913  py1 = py1/pt1;
1914 
1915  TVector *varvec = new TVector(9);
1916  (*varvec)[0] = d0z0[0];
1917  (*varvec)[1] = x1;
1918  (*varvec)[2] = y1;
1919  (*varvec)[3] = px1;
1920  (*varvec)[4] = py1;
1921  (*varvec)[5] = fVtx1->GetX();
1922  (*varvec)[6] = fVtx1->GetY();
1923  (*varvec)[7] = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTPC(trk,AliPID::kProton);
1924  (*varvec)[8] = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTOF(trk,AliPID::kProton);
1925  fProtonCutVarsArray->AddLast(varvec);
1926  }
1927 
1928  Int_t pdgPr = -9999;
1929  if(fUseMCInfo)
1930  {
1931  Int_t labPr = trk->GetLabel();
1932  if(labPr<0) return;
1933  AliAODMCParticle *mcetrk = (AliAODMCParticle*)mcArray->At(labPr);
1934  if(!mcetrk) return;
1935  pdgPr = mcetrk->GetPdgCode();
1936  if(abs(pdgPr)!=2212) return;
1937  fHistoBachPtMCS->Fill(vpr.Pt(),vpr.Rapidity());
1938  }
1939 }
1942 {
1943  //
1945  //
1946 
1947  const char* nameoutput = GetOutputSlot(6)->GetContainer()->GetName();
1948  fV0VariablesTree = new TTree(nameoutput,"v0 variables tree");
1949  Int_t nVar = 3;
1951  TString * fCandidateVariableNames = new TString[nVar];
1952 
1953  fCandidateVariableNames[ 0]="V0Px";
1954  fCandidateVariableNames[ 1]="V0Py";
1955  fCandidateVariableNames[ 2]="V0Pz";
1956 
1957  for (Int_t ivar=0; ivar<nVar; ivar++) {
1958  fV0VariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateV0Variables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
1959  }
1960 
1961  return;
1962 }
1963 //________________________________________________________________________
1964 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillV0ROOTObjects(AliAODv0 *v0, TClonesArray *mcArray)
1965 {
1966  //
1968  //
1969  if(!v0) return;
1970 
1971  TLorentzVector vk0;
1972  vk0.SetXYZM(v0->Px(),v0->Py(),v0->Pz(),v0->MassK0Short());
1973  fHistoK0sMassvsPt->Fill(v0->MassK0Short(),v0->Pt(),vk0.Rapidity());
1974 
1975  if(fDoEventMixing){
1976  Double_t ev0 = sqrt(v0->P2()+pow(v0->MassK0Short(),2));
1977  fV0Tracks->AddLast(new TLorentzVector(v0->Px(),v0->Py(),v0->Pz(),ev0));
1978 
1979  TVector *varvec = new TVector(1);
1980  (*varvec)[0] = v0->DcaV0ToPrimVertex();
1981  fV0CutVarsArray->AddLast(varvec);
1982  }
1983 
1984  Int_t v0pdgcode = -9999;
1985  if(fUseMCInfo)
1986  {
1987  Int_t pdgdgv0[2]={211,211};
1988  Int_t labV0 = v0->MatchToMC(310,mcArray,2,pdgdgv0); // the V0
1989  if(labV0<0) return;
1990  AliAODMCParticle *mcv0trk = (AliAODMCParticle*)mcArray->At(labV0);
1991  if(!mcv0trk) return;
1992  fHistoK0sMassvsPtMCS->Fill(v0->MassK0Short(),v0->Pt(),vk0.Rapidity());
1993  }
1994 }
1995 //_________________________________________________________________
1997  //
1998  // check in which of the pools the current event falls
1999  //
2000 
2001  Int_t theBinZ=TMath::BinarySearch(fNzVtxBins,fZvtxBins,zvert);
2002  if(theBinZ<0 || theBinZ>=fNzVtxBins) return -1;
2003  Int_t theBinM=TMath::BinarySearch(fNCentBins,fCentBins,mult);
2004  if(theBinM<0 || theBinM>=fNCentBins) return -1;
2005  Int_t theBinR=TMath::BinarySearch(fNRPBins,fRPBins,rp);
2006  if(theBinR<0 || theBinR>=fNRPBins) return -1;
2007  return fNRPBins*fNCentBins*theBinZ+fNRPBins*theBinM+theBinR;
2008 }
2009 //_________________________________________________________________
2011  //
2012  // delete the contets of the pool
2013  //
2014  if(poolIndex<0 || poolIndex>=fNOfPools) return;
2015  delete fEventBuffer[poolIndex];
2016  fEventBuffer[poolIndex]=new TTree(Form("EventBuffer_%d",poolIndex), "Temporary buffer for event mixing");
2017 
2018  fEventBuffer[poolIndex]->Branch("zVertex", &fVtxZ);
2019  fEventBuffer[poolIndex]->Branch("centrality", &fCentrality);
2020  fEventBuffer[poolIndex]->Branch("eventInfo", "TObjString",&fEventInfo);
2021  fEventBuffer[poolIndex]->Branch("v1array", "TObjArray", &fV0Tracks);
2022  fEventBuffer[poolIndex]->Branch("v1varsarray", "TObjArray", &fV0CutVarsArray);
2023 
2024  return;
2025 }
2026 //_________________________________________________________________
2028 {
2029  //
2030  // perform mixed event analysis
2031  //
2032 
2033  if(poolIndex<0 || poolIndex>fNzVtxBins*fNCentBins*fNRPBins) return;
2034  if(fEventBuffer[poolIndex]->GetEntries()<fNumberOfEventsForMixing) return;
2035 
2036  Int_t nPro = fProtonTracks->GetEntries();
2037  Int_t nEvents=fEventBuffer[poolIndex]->GetEntries();
2038  Int_t nPro_test = fProtonCutVarsArray->GetEntries();
2039  if(nPro!=nPro_test){
2040  cout<<"Something wrong in mixing machinery"<<endl;
2041  exit(1);
2042  }
2043 
2044  TObjArray* v1array=0x0;
2045  TObjArray* v1varsarray=0x0;
2046  Float_t zVertex,cent;
2047  TObjString* eventInfo=0x0;
2048  fEventBuffer[poolIndex]->SetBranchAddress("eventInfo",&eventInfo);
2049  fEventBuffer[poolIndex]->SetBranchAddress("zVertex", &zVertex);
2050  fEventBuffer[poolIndex]->SetBranchAddress("centrality", &cent);
2051  fEventBuffer[poolIndex]->SetBranchAddress("v1array", &v1array);
2052  fEventBuffer[poolIndex]->SetBranchAddress("v1varsarray", &v1varsarray);
2053  for (Int_t i=0; i<nPro; i++)
2054  {
2055  TLorentzVector* trke=(TLorentzVector*) fProtonTracks->At(i);
2056  if(!trke)continue;
2057  TVector *prvarsarray = (TVector*)fProtonCutVarsArray->At(i);
2058 
2059  for(Int_t iEv=0; iEv<fNumberOfEventsForMixing; iEv++){
2060  fEventBuffer[poolIndex]->GetEvent(iEv + nEvents - fNumberOfEventsForMixing);
2061  Int_t nV01=v1array->GetEntries();
2062  for(Int_t iTr1=0; iTr1<nV01; iTr1++){
2063  TLorentzVector* v01=(TLorentzVector*)v1array->At(iTr1);
2064  if(!v01 ) continue;
2065  if(!fAnalCuts->SelectWithRoughCuts(v01,trke)) continue;
2066 
2067  TVector *v0varsarray = (TVector*) v1varsarray->At(iTr1);
2068  FillMixROOTObjects(trke,v01,prvarsarray,v0varsarray);
2069  }//v0 loop
2070 
2071  }//event loop
2072  }//track loop
2073 }
2074 
2077 {
2081 
2082  const char* nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
2083  fMCVariablesTree = new TTree(nameoutput,"MC variables tree");
2084  Int_t nVar = 1;
2086  TString * fCandidateVariableNames = new TString[nVar];
2087 
2088  fCandidateVariableNames[ 0]="Centrality";
2089 
2090  for (Int_t ivar=0; ivar<nVar; ivar++) {
2091  fMCVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateMCVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
2092  }
2093  return;
2094 }
2096 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillMCROOTObjects(AliAODMCParticle *mcpart, AliAODMCParticle *mcepart, AliAODMCParticle *mcv0part, Int_t decaytype)
2097 {
2098  //
2100  //
2101  if(!mcpart) return;
2102  if(!mcepart) return;
2103  if(!mcv0part) return;
2104 
2105 
2106  if(decaytype==0){
2107  Double_t contmc[3];
2108  contmc[0] = mcpart->Pt();
2109  contmc[1] = mcpart->Y();
2110  contmc[2] = fCentrality;
2111  fHistoLcMCGen->Fill(contmc);
2112  }
2113 }
2116 {
2117  //
2118  // Define proton tree variables
2119  //
2120 
2121  const char* nameoutput = GetOutputSlot(9)->GetContainer()->GetName();
2122  fMCProtonVariablesTree = new TTree(nameoutput,"MC Proton variables tree");
2123  Int_t nVar = 1;
2125  TString * fCandidateVariableNames = new TString[nVar];
2126 
2127  fCandidateVariableNames[ 0]="Centrality";
2128 
2129  for (Int_t ivar=0; ivar<nVar; ivar++) {
2130  fMCProtonVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateMCProtonVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
2131  }
2132  return;
2133 }
2135 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillMCProtonROOTObjects(AliAODMCParticle *mcepart, TClonesArray *mcArray)
2136 {
2137  //
2138  // Fill tree depending on fWriteMCVariableTree
2139  //
2140  if(!mcepart) return;
2141 }
2144 {
2145  //
2146  // Define Mc v0 tree variables
2147  //
2148 
2149  const char* nameoutput = GetOutputSlot(10)->GetContainer()->GetName();
2150  fMCV0VariablesTree = new TTree(nameoutput,"MC v0 variables tree");
2151  Int_t nVar = 1;
2153  TString * fCandidateVariableNames = new TString[nVar];
2154 
2155  fCandidateVariableNames[ 0]="Centrality";
2156 
2157  for (Int_t ivar=0; ivar<nVar; ivar++) {
2158  fMCV0VariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateMCV0Variables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
2159  }
2160  return;
2161 }
2163 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillMCV0ROOTObjects(AliAODMCParticle *mcv0part, TClonesArray *mcArray)
2164 {
2165  //
2166  // Fill histograms or tree depending on fWriteMCVariableTree
2167  //
2168  if(!mcv0part) return;
2169 }
2170 //_________________________________________________________________
2172 {
2173  //
2174  // Analyze AliAODmcparticle
2175  //
2176  Int_t nmcpart = mcArray->GetEntriesFast();
2177 
2178  for(Int_t i=0;i<nmcpart;i++)
2179  {
2180  AliAODMCParticle *mcpart = (AliAODMCParticle*) mcArray->At(i);
2181  if(TMath::Abs(mcpart->GetPdgCode())==4122){
2182  //cout<<"Lambdac"<<endl;
2183  Bool_t p_flag = kFALSE;
2184  Bool_t k0s_flag = kFALSE;
2185  AliAODMCParticle *mcepart = 0;
2186  AliAODMCParticle *mcv0part = 0;
2187  Int_t ndau = mcpart->GetLastDaughter()-mcpart->GetFirstDaughter()+1;
2188  if(ndau==2){
2189  for(Int_t idau=mcpart->GetFirstDaughter();idau<mcpart->GetLastDaughter()+1;idau++)
2190  {
2191  if(idau<0) break;
2192  AliAODMCParticle *mcdau = (AliAODMCParticle*) mcArray->At(idau);
2193  if(!mcdau) continue;
2194  if(TMath::Abs(mcdau->GetPdgCode())==2212){
2195  p_flag = kTRUE;
2196  mcepart = mcdau;
2197  }
2198  if(TMath::Abs(mcdau->GetPdgCode())==311){
2199  k0s_flag = kTRUE;
2200  mcv0part = mcdau;
2201  }
2202  }
2203  }
2204 
2205  Int_t decaytype = -9999;
2206  if(p_flag && k0s_flag) decaytype = 0;
2207 
2208  FillMCROOTObjects(mcpart,mcepart,mcv0part,decaytype);
2209  }
2210 
2211  if(TMath::Abs(mcpart->GetPdgCode())==2212 && mcpart->GetStatus()==1){
2212  AliESDtrackCuts *esdcuts = fAnalCuts->GetTrackCuts();
2214  esdcuts->GetEtaRange(etamin,etamax);
2215  if(fabs(mcpart->Eta())<etamax){
2216  fHistoBachPtMCGen->Fill(mcpart->Pt(),mcpart->Y());
2217  }
2218  FillMCProtonROOTObjects(mcpart, mcArray);
2219  }
2220  if(TMath::Abs(mcpart->GetPdgCode())==321 && mcpart->GetStatus()==1){
2221  AliESDtrackCuts *esdcuts = fAnalCuts->GetTrackCuts();
2223  esdcuts->GetEtaRange(etamin,etamax);
2224  if(fabs(mcpart->Eta())<etamax){
2225  fHistoKaonPtMCGen->Fill(mcpart->Pt(),mcpart->Y());
2226  }
2227  }
2228  if(TMath::Abs(mcpart->GetPdgCode())==310){
2229  Double_t etamin, etamax, rapmin, rapmax;
2230  fAnalCuts->GetProdV0EtaRange(etamin,etamax);
2231  fAnalCuts->GetProdV0RapRange(rapmin,rapmax);
2232  if((fabs(mcpart->Y())<rapmax) && (fabs(mcpart->Eta())<etamax)){
2233  fHistoK0sMassvsPtMCGen->Fill(0.497, mcpart->Pt(),mcpart->Y());
2234  }
2235  FillMCV0ROOTObjects(mcpart, mcArray);
2236  }
2237  }
2238 
2239  return kTRUE;
2240 }
2241 //________________________________________________________________________
2242 Int_t AliAnalysisTaskSELc2pK0sfromAODtracks::MatchToMC(AliAODRecoCascadeHF *elobj, TClonesArray *mcArray, Int_t *pdgarray_pr, Int_t *pdgarray_v0, Int_t *labelarray_pr, Int_t *labelarray_v0, Int_t &ngen_pr, Int_t &ngen_v0)
2243 {
2244  //
2245  // Match to MC
2246  //
2247  for(Int_t i=0;i<100;i++){
2248  pdgarray_pr[i] = -9999;
2249  labelarray_pr[i] = -9999;
2250  pdgarray_v0[i] = -9999;
2251  labelarray_v0[i] = -9999;
2252  }
2253  ngen_pr = 0;
2254  ngen_v0 = 0;
2255 
2256  AliVTrack *trk = dynamic_cast<AliVTrack*>(elobj->GetBachelor());
2257  if(!trk) return -1;
2258  Int_t labPr = trk->GetLabel();
2259  if(labPr<0) return -1;
2260  AliAODMCParticle *mcetrk = (AliAODMCParticle*)mcArray->At(labPr);
2261  if(!mcetrk) return -1;
2262  labelarray_pr[0] = labPr;
2263  pdgarray_pr[0] = mcetrk->GetPdgCode();
2264  ngen_pr ++;
2265 
2266  AliAODMCParticle *mcprimpr=0;
2267  mcprimpr = mcetrk;
2268  while(mcprimpr->GetMother()>=0) {
2269  Int_t labprim_pr=mcprimpr->GetMother();
2270  AliAODMCParticle *tmcprimpr = (AliAODMCParticle*)mcArray->At(labprim_pr);
2271  if(!tmcprimpr) {
2272  break;
2273  }
2274 
2275  mcprimpr = tmcprimpr;
2276  pdgarray_pr[ngen_pr] = mcprimpr->GetPdgCode();
2277  labelarray_pr[ngen_pr] = labprim_pr;
2278  ngen_pr ++;
2279  if(ngen_pr==100) break;
2280  }
2281 
2282  AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(elobj->Getv0());
2283  if(!theV0) return -1;
2284  Int_t pdgdgv0[2]={211,211};
2285  Int_t labV0 = theV0->MatchToMC(310,mcArray,2,pdgdgv0); // the V0
2286  if(labV0<0) return -1;
2287  AliAODMCParticle *mcv0 = (AliAODMCParticle*)mcArray->At(labV0);
2288  if(!mcv0) return -1;
2289  labelarray_v0[0] = labV0;
2290  pdgarray_v0[0] = mcv0->GetPdgCode();
2291  ngen_v0 ++;
2292 
2293  AliAODMCParticle *mcprimv0=0;
2294  mcprimv0 = mcv0;
2295  while(mcprimv0->GetMother()>=0) {
2296  Int_t labprim_v0=mcprimv0->GetMother();
2297  AliAODMCParticle *tmcprimv0 = (AliAODMCParticle*)mcArray->At(labprim_v0);
2298  if(!tmcprimv0) {
2299  break;
2300  }
2301 
2302  mcprimv0 = tmcprimv0;
2303  pdgarray_v0[ngen_v0] = mcprimv0->GetPdgCode();
2304  labelarray_v0[ngen_v0] = labprim_v0;
2305  ngen_v0 ++;
2306  if(ngen_v0==100) break;
2307  }
2308 
2309  Bool_t same_flag = kFALSE;
2310  Int_t matchedlabel=-9999;
2311  for(Int_t iemc=0;iemc<ngen_pr;iemc++){
2312  for(Int_t ivmc=0;ivmc<ngen_v0;ivmc++){
2313  if(labelarray_pr[iemc]==labelarray_v0[ivmc]){
2314  same_flag = kTRUE;
2315  matchedlabel = labelarray_pr[iemc];
2316  break;
2317  }
2318  }
2319  if(same_flag) break;
2320  }
2321 
2322  return matchedlabel;
2323 
2324 }
2325 
2326 //________________________________________________________________________
2328  //
2329  // Stores the pointer to the global track
2330  // copied from femtoscopy/k0analysis/plamanalysis
2331  //
2332 
2333  // Check that the id is positive
2334  if(track->GetID()<0){
2335  // printf("Warning: track has negative ID: %d\n",track->GetID());
2336  return;
2337  }
2338 
2339  // Check id is not too big for buffer
2340  if(track->GetID()>=fTrackBuffSize){
2341  printf("Warning: track ID too big for buffer: ID: %d, buffer %d\n"
2342  ,track->GetID(),fTrackBuffSize);
2343  return;
2344  }
2345 
2346  // Warn if we overwrite a track
2347  if(fGTI[track->GetID()]){
2348  // Seems like there are FilterMap 0 tracks
2349  // that have zero TPCNcls, don't store these!
2350  if( (!track->GetFilterMap()) &&
2351  (!track->GetTPCNcls()) )
2352  return;
2353 
2354  // Imagine the other way around,
2355  // the zero map zero clusters track
2356  // is stored and the good one wants
2357  // to be added. We ommit the warning
2358  // and just overwrite the 'bad' track
2359  if( fGTI[track->GetID()]->GetFilterMap() ||
2360  fGTI[track->GetID()]->GetTPCNcls() ){
2361  // If we come here, there's a problem
2362  printf("Warning! global track info already there!");
2363  printf(" TPCNcls track1 %u track2 %u",
2364  (fGTI[track->GetID()])->GetTPCNcls(),track->GetTPCNcls());
2365  printf(" FilterMap track1 %u track2 %u\n",
2366  (fGTI[track->GetID()])->GetFilterMap(),track->GetFilterMap());
2367  }
2368  } // Two tracks same id
2369 
2370  // // There are tracks with filter bit 0,
2371  // // do they have TPCNcls stored?
2372  // if(!track->GetFilterMap()){
2373  // printf("Filter map is zero, TPCNcls: %u\n"
2374  // ,track->GetTPCNcls());
2375  // }
2376 
2377  // Assign the pointer
2378  (fGTI[track->GetID()]) = track;
2379  (fGTIndex[track->GetID()]) = index;
2380 }
2381 //________________________________________________________________________
2383  // Sets all the pointers to zero. To be called at
2384  // the beginning or end of an event
2385  for(UShort_t i=0;i<fTrackBuffSize;i++){
2386  fGTI[i]=0;
2387  fGTIndex[i]=-9999;
2388  }
2389 }
Int_t charge
TTree * fMCVariablesTree
! tree of the candidate variables after track selection on output slot 7
void FillProtonROOTObjects(AliAODTrack *trk, TClonesArray *mcArray)
Bool_t SingleTrkCuts(AliAODTrack *trk, AliAODTrack *trkpid, AliAODVertex *vtx)
double Double_t
Definition: External.C:58
Definition: External.C:260
Int_t MatchToMC(AliAODRecoCascadeHF *elobj, TClonesArray *mcArray, Int_t *pdgarray_pr, Int_t *pdgarray_v0, Int_t *labelarray_pr, Int_t *labelarray_v0, Int_t &ngen_pr, Int_t &ngen_v0)
Definition: External.C:236
TH2F * fHistoBachPtMCGen
! Bachelor pT histogram (efficiency denominator)
TObjArray * fProtonTracks
unique event id for mixed event check
Bool_t GetUseCombined()
Definition: AliAODPidHF.h:165
void SelectTrack(const AliVEvent *event, Int_t trkEntries, Int_t &nSeleTrks, Int_t *seleFlags, TClonesArray *mcArray)
AliAODVertex * PrimaryVertex(const TObjArray *trkArray, AliVEvent *event)
TH1F * fHistoProbProton
! Probability to be proton histogram
AliAODv0 * Getv0() const
TH1F * fHReactionPlane
! Histogram to check Reaction plane
char Char_t
Definition: External.C:18
Float_t * fCandidateProtonVariables
! variables to be written to the tree
TTree * fVariablesTree
flag to decide whether to write the candidate variables on a tree variables
void GetProdV0EtaRange(Double_t &a, Double_t &b)
TH2F * fHistoKaonPtMCGen
! Kaon pT histogram (efficiency denominator)
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:256
TH2F * fHistoKaonPtMCS
! Kaon pT histogram (efficiency numerator)
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
Bool_t SingleV0Cuts(AliAODv0 *v0, AliAODVertex *vert)
Double_t InvMassLctoK0sP() const
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:246
Float_t * fCandidateV0Variables
! variables to be written to the tree
const Double_t etamin
TH1F * fHCentrality
! Histogram to check Centrality
TTree * fMCProtonVariablesTree
! tree of the candidate variables after track selection on output slot 9
int Int_t
Definition: External.C:63
TTree * fMCV0VariablesTree
! tree of the candidate variables after track selection on output slot 10
const UShort_t fTrackBuffSize
Array of integers to keep the index of tpc only track.
float Float_t
Definition: External.C:68
AliAODRecoCascadeHF * MakeCascadeHF(AliAODv0 *casc, AliAODTrack *trk, AliAODTrack *trkpid, AliAODEvent *aod, AliAODVertex *vert)
TTree * fProtonVariablesTree
! tree of the candidate variables after track selection on output slot 5
Double_t CalculateLcCosPAXY(AliAODRecoDecayHF *obj)
virtual void UserCreateOutputObjects()
Implementation of interface methods.
void MakeAnalysis(AliAODEvent *aod, TClonesArray *mcArray)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Definition: External.C:228
AliAODTrack * GetBachelor() const
TH2F * fHistoBachPtMCS
! Bachelor pT histogram (efficiency numerator)
TObjArray * fV0CutVarsArray
array of RDHF cut information
AliESDtrackCuts * GetTrackCuts() const
Definition: AliRDHFCuts.h:261
void FillMCV0ROOTObjects(AliAODMCParticle *mcv0part, TClonesArray *mcArray)
void SetProngIDs(Int_t nIDs, UShort_t *id)
void GetProdV0RapRange(Double_t &a, Double_t &b)
Double_t nEvents
plot quality messages
Bool_t fWriteMCVariableTree
flag to decide whether to write the candidate variables on a tree variables
AliPIDResponse * GetPidResponse() const
Definition: AliAODPidHF.h:160
Bool_t fWriteEachVariableTree
flag to decide whether to write the candidate variables on a tree variables
Int_t GetPoolIndex(Double_t zvert, Double_t mult, Double_t rp)
void SelectV0(const AliVEvent *event, Int_t nV0, Int_t &nSeleV0, Bool_t *seleV0Flags, TClonesArray *mcArray)
Double_t CosV0PointingAngle() 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)
Bool_t IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:345
short Short_t
Definition: External.C:23
void FillROOTObjects(AliAODRecoCascadeHF *lcobj, AliAODv0 *v0, AliAODTrack *trk, AliAODTrack *trkpid, AliAODEvent *aod, TClonesArray *mcarray)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Double_t DecayLengthXY() const
TTree * fV0VariablesTree
! tree of the candidate variables after track selection on output slot 6
Int_t fNzVtxBins
maximum number of events to be used in event mixing
void FillV0ROOTObjects(AliAODv0 *v0, TClonesArray *mcArray)
Float_t * fCandidateMCVariables
! variables to be written to the tree
void FillMixROOTObjects(TLorentzVector *pt, TLorentzVector *ev, TVector *tinfo, TVector *vinfo)
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
void FillMCROOTObjects(AliAODMCParticle *part, AliAODMCParticle *mcepart, AliAODMCParticle *mcv0part, Int_t decaytype)
Bool_t fIsMB
Reconstruct primary vertex excluding candidate tracks.
const Int_t nVar
TList * fOutputAll
! User Output slot 3 //analysis histograms
Bool_t TagV0(AliAODTrack *etrk, AliAODEvent *evt, Int_t ntrk, Double_t &minmass)
const Double_t etamax
Float_t * fCandidateMCV0Variables
! variables to be written to the tree
Bool_t SelectWithRoughCuts(AliAODv0 *v0, AliAODTrack *trk1)
unsigned short UShort_t
Definition: External.C:28
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:269
TH1F * fCEvents
! Histogram to check selected events
const char Option_t
Definition: External.C:48
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:315
Float_t * fCandidateMCProtonVariables
! variables to be written to the tree
void FillMCProtonROOTObjects(AliAODMCParticle *mcepart, TClonesArray *mcArray)
bool Bool_t
Definition: External.C:53
Double_t fCentBins[100]
number of centrality bins
AliNormalizationCounter * fCounter
EvNumber counter.
AliAODVertex * CallPrimaryVertex(AliAODv0 *v0, AliAODTrack *trk, AliAODEvent *evt)
TObjArray * fV0Tracks
array of electron-compatible tracks
Int_t * fGTIndex
Array of pointers, just nicely sorted according to the id.
AliAODVertex * ReconstructSecondaryVertex(AliAODv0 *casc, AliAODTrack *trk, AliAODEvent *aod)
Int_t fDoEventMixing
Size of the above array, ~12000 for PbPb.
TH1F * fHistoV0CosPA
! V0 cosine pointing angle to primary vertex
TObjArray * fProtonCutVarsArray
array of lambda-compatible tracks