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