AliPhysics  vAN-20150924 (e816f45)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 <THnSparse.h>
45 #include <TLorentzVector.h>
46 #include <TTree.h>
47 #include "TROOT.h"
48 #include <TDatabasePDG.h>
49 #include <AliAnalysisDataSlot.h>
50 #include <AliAnalysisDataContainer.h>
51 #include "AliStack.h"
52 #include "AliMCEvent.h"
53 #include "AliAnalysisManager.h"
54 #include "AliAODMCHeader.h"
55 #include "AliAODHandler.h"
56 #include "AliLog.h"
57 #include "AliExternalTrackParam.h"
58 #include "AliAODVertex.h"
59 #include "AliESDVertex.h"
60 #include "AliAODRecoDecay.h"
61 #include "AliAODRecoDecayHF.h"
62 #include "AliAODRecoCascadeHF.h"
63 #include "AliESDtrack.h"
64 #include "AliAODTrack.h"
65 #include "AliAODv0.h"
66 #include "AliAODcascade.h"
67 #include "AliAODMCParticle.h"
68 #include "AliAnalysisTaskSE.h"
70 #include "AliPIDResponse.h"
71 #include "AliPIDCombined.h"
72 #include "AliTOFPIDResponse.h"
73 #include "AliAODPidHF.h"
74 #include "AliInputEventHandler.h"
75 #include "AliESDtrackCuts.h"
76 #include "AliNeutralTrackParam.h"
77 #include "AliKFParticle.h"
78 #include "AliKFVertex.h"
79 #include "AliExternalTrackParam.h"
80 #include "AliESDtrack.h"
81 #include "AliCentrality.h"
82 #include "AliVertexerTracks.h"
83 
84 using std::cout;
85 using std::endl;
86 
90 
91 //__________________________________________________________________________
93  AliAnalysisTaskSE(),
94  fUseMCInfo(kFALSE),
95  fOutput(0),
96  fOutputAll(0),
97  fListCuts(0),
98  fCEvents(0),
99  fHTrigger(0),
100  fHCentrality(0),
101  fAnalCuts(0),
102  fIsEventSelected(kFALSE),
103  fWriteVariableTree(kFALSE),
104  fVariablesTree(0),
105  fReconstructPrimVert(kFALSE),
106  fIsMB(kFALSE),
107  fIsSemi(kFALSE),
108  fIsCent(kFALSE),
109  fIsINT7(kFALSE),
110  fIsEMC7(kFALSE),
111  fCandidateVariables(),
112  fVtx1(0),
113  fV1(0),
114  fBzkG(0),
115  fCentrality(0),
116  fTriggerCheck(0),
117  fHistoLcK0SpMass(0),
118  fHistoBachPt(0),
119  fHistod0Bach(0),
120  fHistod0V0(0),
121  fHistod0d0(0),
122  fHistoV0CosPA(0),
123  fHistoProbProton(0),
124  fHistoDecayLength(0),
125  fHistoK0SMass(0)
126 {
127  //
128  // Default Constructor.
129  //
130 }
131 
132 //___________________________________________________________________________
135  Bool_t writeVariableTree) :
136  AliAnalysisTaskSE(name),
137  fUseMCInfo(kFALSE),
138  fOutput(0),
139  fOutputAll(0),
140  fListCuts(0),
141  fCEvents(0),
142  fHTrigger(0),
143  fHCentrality(0),
144  fAnalCuts(analCuts),
145  fIsEventSelected(kFALSE),
146  fWriteVariableTree(writeVariableTree),
147  fVariablesTree(0),
148  fReconstructPrimVert(kFALSE),
149  fIsMB(kFALSE),
150  fIsSemi(kFALSE),
151  fIsCent(kFALSE),
152  fIsINT7(kFALSE),
153  fIsEMC7(kFALSE),
154  fCandidateVariables(),
155  fVtx1(0),
156  fV1(0),
157  fBzkG(0),
158  fCentrality(0),
159  fTriggerCheck(0),
160  fHistoLcK0SpMass(0),
161  fHistoBachPt(0),
162  fHistod0Bach(0),
163  fHistod0V0(0),
164  fHistod0d0(0),
165  fHistoV0CosPA(0),
166  fHistoProbProton(0),
167  fHistoDecayLength(0),
168  fHistoK0SMass(0)
169 {
170  //
171  // Constructor. Initialization of Inputs and Outputs
172  //
173  Info("AliAnalysisTaskSELc2pK0sfromAODtracks","Calling Constructor");
174 
175  DefineOutput(1,TList::Class()); //conters
176  DefineOutput(2,TList::Class());
177  if(writeVariableTree){
178  DefineOutput(3,TTree::Class()); //My private output
179  }else{
180  DefineOutput(3,TList::Class()); //conters
181  }
182 }
183 
184 //___________________________________________________________________________
186  //
187  // destructor
188  //
189  Info("~AliAnalysisTaskSELc2pK0sfromAODtracks","Calling Destructor");
190 
191  if (fOutput) {
192  delete fOutput;
193  fOutput = 0;
194  }
195 
196  if (fOutputAll) {
197  delete fOutputAll;
198  fOutputAll = 0;
199  }
200 
201  if (fListCuts) {
202  delete fListCuts;
203  fListCuts = 0;
204  }
205 
206 
207  if (fAnalCuts) {
208  delete fAnalCuts;
209  fAnalCuts = 0;
210  }
211 
212  if (fVariablesTree) {
213  delete fVariablesTree;
214  fVariablesTree = 0;
215  }
216 
217 }
218 
219 //_________________________________________________
221  //
222  // Initialization
223  //
224  //
225 
226  fIsEventSelected=kFALSE;
227 
228  if (fDebug > 1) AliInfo("Init");
229 
230  fListCuts = new TList();
231  fListCuts->SetOwner();
232  fListCuts->SetName("ListCuts");
234  PostData(2,fListCuts);
235 
236  return;
237 }
238 
239 //_________________________________________________
241 {
242  //
243  // UserExec
244  //
245 
246  if (!fInputEvent) {
247  AliError("NO EVENT FOUND!");
248  return;
249  }
250  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
251 
252  fCEvents->Fill(1);
253 
254  //------------------------------------------------
255  // First check if the event has proper vertex and B
256  //------------------------------------------------
257 
258  fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
259  if (!fVtx1) return;
260 
261  Double_t pos[3],cov[6];
262  fVtx1->GetXYZ(pos);
263  fVtx1->GetCovarianceMatrix(cov);
264  fV1 = new AliESDVertex(pos,cov,100.,100,fVtx1->GetName());
265 
266  fBzkG = (Double_t)aodEvent->GetMagneticField();
267  AliKFParticle::SetField(fBzkG);
268  if (TMath::Abs(fBzkG)<0.001) {
269  delete fV1;
270  return;
271  }
272  fCEvents->Fill(2);
273 
274  //------------------------------------------------
275  // Event selection
276  //------------------------------------------------
277  Bool_t fIsTriggerNotOK = fAnalCuts->IsEventRejectedDueToTrigger();
278  if(!fIsTriggerNotOK) fCEvents->Fill(3);
279  fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent); // better to initialize before CheckEventSelection call
280  if(!fIsEventSelected) {
281  delete fV1;
282  return;
283  }
284  fCEvents->Fill(4);
285 
286  fIsMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB)==(AliVEvent::kMB);
287  fIsSemi=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kSemiCentral)==(AliVEvent::kSemiCentral);
288  fIsCent=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kCentral)==(AliVEvent::kCentral);
289  fIsINT7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kINT7)==(AliVEvent::kINT7);
290  fIsEMC7=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kEMC7)==(AliVEvent::kEMC7);
292  if(fIsMB) fHTrigger->Fill(1);
293  if(fIsSemi) fHTrigger->Fill(2);
294  if(fIsCent) fHTrigger->Fill(3);
295  if(fIsINT7) fHTrigger->Fill(4);
296  if(fIsEMC7) fHTrigger->Fill(5);
297  if(fIsMB|fIsSemi|fIsCent) fHTrigger->Fill(7);
298  if(fIsINT7|fIsEMC7) fHTrigger->Fill(8);
299  if(fIsMB&fIsSemi) fHTrigger->Fill(10);
300  if(fIsMB&fIsCent) fHTrigger->Fill(11);
301  if(fIsINT7&fIsEMC7) fHTrigger->Fill(12);
302 
303  AliCentrality *cent = aodEvent->GetCentrality();
304  fCentrality = cent->GetCentralityPercentile("V0M");
305  fHCentrality->Fill(fCentrality);
306 
307  //------------------------------------------------
308  // Check if the event has v0 candidate
309  //------------------------------------------------
310  //Int_t nv0 = aodEvent->GetNumberOfV0s();
311  fCEvents->Fill(5);
312 
313  //------------------------------------------------
314  // MC analysis setting
315  //------------------------------------------------
316  TClonesArray *mcArray = 0;
317  AliAODMCHeader *mcHeader=0;
318  if (fUseMCInfo) {
319  // MC array need for maching
320  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
321  if (!mcArray) {
322  AliError("Could not find Monte-Carlo in AOD");
323  return;
324  }
325  fCEvents->Fill(6); // in case of MC events
326 
327  // load MC header
328  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
329  if (!mcHeader) {
330  AliError("AliAnalysisTaskSELc2pK0sfromAODtracks::UserExec: MC header branch not found!\n");
331  return;
332  }
333  fCEvents->Fill(7); // in case of MC events
334 
335  Double_t zMCVertex = mcHeader->GetVtxZ();
336  if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()) {
337  AliDebug(2,Form("Event rejected: abs(zVtxMC)=%f > fAnalCuts->GetMaxVtxZ()=%f",zMCVertex,fAnalCuts->GetMaxVtxZ()));
338  return;
339  } else {
340  fCEvents->Fill(17); // in case of MC events
341  }
342  }
343 
344  //------------------------------------------------
345  // Main analysis done in this function
346  //------------------------------------------------
347  MakeAnalysis(aodEvent,mcArray);
348 
349 
350  PostData(1,fOutput);
351  if(fWriteVariableTree){
352  PostData(3,fVariablesTree);
353  }else{
354  PostData(3,fOutputAll);
355  }
356 
357  fIsEventSelected=kFALSE;
358 
359  delete fV1;
360  return;
361 }
362 
363 //________________________________________ terminate ___________________________
365 {
366  // The Terminate() function is the last function to be called during
367  // a query. It always runs on the client, it can be used to present
368  // the results graphically or save the results to file.
369 
370  //AliInfo("Terminate","");
371  AliAnalysisTaskSE::Terminate();
372 
373  fOutput = dynamic_cast<TList*> (GetOutputData(1));
374  if (!fOutput) {
375  AliError("fOutput not available");
376  return;
377  }
378 
379  if(!fWriteVariableTree){
380  fOutputAll = dynamic_cast<TList*> (GetOutputData(3));
381  if (!fOutputAll) {
382  AliError("fOutputAll not available");
383  return;
384  }
385  }
386 
387  return;
388 }
389 
390 //___________________________________________________________________________
392 {
393  //
394  // UserCreateOutputObject
395  //
396  //AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
397 
398  //------------------------------------------------
399  // output object setting
400  //------------------------------------------------
401  fOutput = new TList();
402  fOutput->SetOwner();
403  fOutput->SetName("chist0");
404  DefineGeneralHistograms(); // define general histograms
405  PostData(1,fOutput);
406 
408  if (fWriteVariableTree) {
409  PostData(3,fVariablesTree);
410  }else{
411  fOutputAll = new TList();
412  fOutputAll->SetOwner();
413  fOutputAll->SetName("anahisto");
414  DefineAnalysisHistograms(); // define general histograms
415  PostData(3,fOutputAll);
416  }
417 
418 
419  return;
420 }
421 
422 //-------------------------------------------------------------------------------
424 (
425  AliAODEvent *aodEvent, TClonesArray *mcArray
426  )
427 {
428  //
429  // Main Analysis part
430  //
431 
432  Int_t nV0s= aodEvent->GetNumberOfV0s();
433  if (nV0s==0) {
434  return;
435  }
436  Int_t nTracks= aodEvent->GetNumberOfTracks();
437  if (nTracks==0) {
438  return;
439  }
440 
441  //------------------------------------------------
442  // Arrays to store MC matching information
443  //------------------------------------------------
444  Int_t usedmclab[20];//Used Lc Label: Assuming there are less than 20 Lc/evt
445  Int_t nusedmclab[20];//Number of times the Lc label is used: Assuming there are less than 20 Lc/evt
446  for(Int_t i=0;i<20;i++) {
447  usedmclab[i]=-9999;
448  nusedmclab[i]=0;
449  }
450 
451  //------------------------------------------------
452  // V0 loop
453  //------------------------------------------------
454  for (Int_t iv0 = 0; iv0<nV0s; iv0++) {
455  AliAODv0 *v0 = aodEvent->GetV0(iv0);
456  if(!v0) continue;
457  if(!fAnalCuts->SingleV0Cuts(v0,fVtx1)) continue;
458 
459  AliAODTrack *cptrack = (AliAODTrack*)(v0->GetDaughter(0));
460  AliAODTrack *cntrack = (AliAODTrack*)(v0->GetDaughter(1));
461 
462  //------------------------------------------------
463  // track loop
464  //------------------------------------------------
465  for (Int_t itrk = 0; itrk<nTracks; itrk++) {
466  AliAODTrack *trk = (AliAODTrack*)aodEvent->GetTrack(itrk);
467  if(trk->GetID()<0) continue;
468  if(!fAnalCuts->SingleTrkCuts(trk)) continue;
469 
470  Int_t cpid = cptrack->GetID();
471  Int_t cnid = cntrack->GetID();
472  Int_t lpid = trk->GetID();
473  if((cpid==lpid)||(cnid==lpid)) continue;
474 
475  if(!fAnalCuts->SelectWithRoughCuts(v0,trk)) continue;
476 
477  AliAODVertex *secVert = ReconstructSecondaryVertex(v0,trk,aodEvent);
478  if(!secVert) continue;
479 
480  AliAODRecoCascadeHF *lcobj = MakeCascadeHF(v0,trk,aodEvent,secVert);
481  if(!lcobj) {
482  continue;
483  }
484 
485  AliAODMCParticle *mclc = 0;
486  AliAODMCParticle *mcproton = 0;
487  AliAODMCParticle *mck0s = 0;
488  Int_t mclablc = 0;
489  Int_t nmclablc = 0;
490  if(fUseMCInfo)
491  {
492  Int_t pdgDg[2]={2212,310};
493  Int_t pdgDgv0[2]={211,211};
494  mclablc = lcobj->MatchToMC(4122,pdgDg[1],pdgDg,pdgDgv0,mcArray,kTRUE);
495  if(mclablc>-1){
496  mclc = (AliAODMCParticle*) mcArray->At(mclablc);
497  for(Int_t ia=0;ia<20;ia++){
498  if(usedmclab[ia]==mclablc){
499  nusedmclab[ia]++;
500  nmclablc = nusedmclab[ia];
501  break;
502  }
503  if(usedmclab[ia]==-9999){
504  usedmclab[ia]=mclablc;
505  nusedmclab[ia]++;
506  nmclablc = nusedmclab[ia];
507  break;
508  }
509  }
510  Int_t mcprotonlabel = mclc->GetDaughter(0);
511  if(mcprotonlabel>=0)
512  mcproton=(AliAODMCParticle*) mcArray->At(mcprotonlabel);
513  Int_t mck0slabel = mclc->GetDaughter(1);
514  if(mck0slabel>=0)
515  mck0s=(AliAODMCParticle*) mcArray->At(mck0slabel);
516  }
517  }
518 
519  FillROOTObjects(lcobj,mclc,mcproton,mck0s,nmclablc);
520 
521  lcobj->GetSecondaryVtx()->RemoveDaughters();
522  lcobj->UnsetOwnPrimaryVtx();
523  delete lcobj;lcobj=NULL;
524  delete secVert;
525  }
526  }
527 
528 }
529 
531 void AliAnalysisTaskSELc2pK0sfromAODtracks::FillROOTObjects(AliAODRecoCascadeHF *lcobj, AliAODMCParticle *mcpart, AliAODMCParticle *mcproton, AliAODMCParticle *mck0s, Int_t mcnused)
532 {
533  //
534  // Fill histograms or tree depending on fWriteVariableTree
535  //
536 
537  AliAODTrack *trk = lcobj->GetBachelor();
538  AliAODv0 *v0 = lcobj->Getv0();
539 
540  fCandidateVariables[ 0] = lcobj->InvMassLctoK0sP();
541  fCandidateVariables[ 1] = lcobj->Px();
542  fCandidateVariables[ 2] = lcobj->Py();
543  fCandidateVariables[ 3] = lcobj->Pz();
544  fCandidateVariables[ 4] = v0->MassK0Short();
545  fCandidateVariables[ 5] = lcobj->PxProng(0);
546  fCandidateVariables[ 6] = lcobj->PyProng(0);
547  fCandidateVariables[ 7] = lcobj->PzProng(0);
548  fCandidateVariables[ 8] = lcobj->PxProng(1);
549  fCandidateVariables[ 9] = lcobj->PyProng(1);
550  fCandidateVariables[10] = lcobj->PzProng(1);
551  fCandidateVariables[11] = fVtx1->GetX();
552  fCandidateVariables[12] = fVtx1->GetY();
553  fCandidateVariables[13] = fVtx1->GetZ();
555  fCandidateVariables[15] = lcobj->DecayLengthXY();
556  fCandidateVariables[16] = (Float_t) fAnalCuts->CalculateLcCosPAXY(lcobj);
557 
558  Double_t nSigmaTPCpr=-9999.;
559  Double_t nSigmaTOFpr=-9999.;
560  Double_t probProton=-9999.;
561  if(fAnalCuts->GetIsUsePID())
562  {
563  nSigmaTPCpr = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTPC(trk,AliPID::kProton);
564  nSigmaTOFpr = fAnalCuts->GetPidHF()->GetPidResponse()->NumberOfSigmasTOF(trk,AliPID::kProton);
566  probProton = fAnalCuts->GetProtonProbabilityTPCTOF(trk);
567  }
568  fCandidateVariables[17] = nSigmaTPCpr;
569  fCandidateVariables[18] = nSigmaTOFpr;
570  fCandidateVariables[19] = probProton;
571  }
572 
573  fCandidateVariables[20] = -9999;
574  fCandidateVariables[21] = -9999;
575  fCandidateVariables[22] = -9999;
576  fCandidateVariables[23] = -9999;
577  fCandidateVariables[24] = -9999;
578  fCandidateVariables[25] = -9999;
579  fCandidateVariables[26] = -9999;
580  if(fUseMCInfo){
581  if(mcpart){
582  fCandidateVariables[20] = mcpart->Label();
583  fCandidateVariables[21] = mcnused;
584  fCandidateVariables[22] = mcpart->GetPdgCode();
585  Double_t mcprimvertx = mcpart->Xv();
586  Double_t mcprimverty = mcpart->Yv();
587  Double_t mcsecvertx = mcproton->Xv();
588  Double_t mcsecverty = mcproton->Yv();
589  Double_t recosecvertx = lcobj->GetSecondaryVtx()->GetX();
590  Double_t recosecverty = lcobj->GetSecondaryVtx()->GetY();
591  fCandidateVariables[23] = TMath::Sqrt((mcsecvertx-mcprimvertx)*(mcsecvertx-mcprimvertx)+(mcsecverty-mcprimverty)*(mcsecverty-mcprimverty));
592  fCandidateVariables[24] = TMath::Sqrt((recosecvertx-mcprimvertx)*(recosecvertx-mcprimvertx)+(recosecverty-mcprimverty)*(recosecverty-mcprimverty));
593  Double_t vecx_vert = recosecvertx-mcprimvertx;
594  Double_t vecy_vert = recosecverty-mcprimverty;
595  Double_t vecl_vert = TMath::Sqrt(vecx_vert*vecx_vert+vecy_vert*vecy_vert);
596  Double_t vecx_mom = lcobj->Px();
597  Double_t vecy_mom = lcobj->Py();
598  Double_t vecl_mom = lcobj->Pt();
599  if(vecl_vert>0.&&vecl_mom>0.)
600  fCandidateVariables[25] = (vecx_vert*vecx_mom+vecy_vert*vecy_mom)/vecl_vert/vecl_mom;
601  fCandidateVariables[26] = mcpart->Pt();
602  }
603  }
604 
606  fVariablesTree->Fill();
607  else{
609  {
610  Double_t cont[3];
611  cont[0] = lcobj->InvMassLctoK0sP();
612  cont[1] = lcobj->Pt();
613  cont[2] = fCentrality;
614  fHistoLcK0SpMass->Fill(cont);
615 
616  fHistoBachPt->Fill(trk->Pt());
617  fHistod0Bach->Fill(lcobj->Getd0Prong(0));
618  fHistod0V0->Fill(lcobj->Getd0Prong(1));
619  fHistod0d0->Fill(lcobj->Getd0Prong(0)*lcobj->Getd0Prong(1));
620  fHistoV0CosPA->Fill(lcobj->CosV0PointingAngle());
621  fHistoProbProton->Fill(probProton);
623  fHistoK0SMass->Fill(v0->MassK0Short());
624  }
625  }
626  return;
627 }
628 
631 {
632  //
633  // Define tree variables
634  //
635 
636  const char* nameoutput = GetOutputSlot(3)->GetContainer()->GetName();
637  fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
638  Int_t nVar = 27;
639  fCandidateVariables = new Float_t [nVar];
640  TString * fCandidateVariableNames = new TString[nVar];
641 
642  fCandidateVariableNames[ 0]="InvMassLc2pK0s";
643  fCandidateVariableNames[ 1]="LcPx";
644  fCandidateVariableNames[ 2]="LcPy";
645  fCandidateVariableNames[ 3]="LcPz";
646  fCandidateVariableNames[ 4]="massK0Short";
647  fCandidateVariableNames[ 5]="V0Px";
648  fCandidateVariableNames[ 6]="V0Py";
649  fCandidateVariableNames[ 7]="V0Pz";
650  fCandidateVariableNames[ 8]="BachPx";
651  fCandidateVariableNames[ 9]="BachPy";
652  fCandidateVariableNames[10]="BachPz";
653  fCandidateVariableNames[11]="PrimVertx";
654  fCandidateVariableNames[12]="PrimVerty";
655  fCandidateVariableNames[13]="PrimVertz";
656  fCandidateVariableNames[14]="Centrality";
657  fCandidateVariableNames[15]="DecayLengthXY";
658  fCandidateVariableNames[16]="LcCosPAXY";
659  fCandidateVariableNames[17]="nSigmaTPCpr";
660  fCandidateVariableNames[18]="nSigmaTOFpr";
661  fCandidateVariableNames[19]="probProton";
662  fCandidateVariableNames[20]="mclcID";
663  fCandidateVariableNames[21]="mcnused";
664  fCandidateVariableNames[22]="mcpdgcode";
665  fCandidateVariableNames[23]="mcdecaylength";
666  fCandidateVariableNames[24]="mcdecaylength_secsmear";
667  fCandidateVariableNames[25]="mclccospaxy";
668  fCandidateVariableNames[26]="mclcpt";
669 
670  for (Int_t ivar=0; ivar<nVar; ivar++) {
671  fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
672  }
673 
674  return;
675 }
676 
679  //
680  // This is to define general histograms
681  //
682 
683  fCEvents = new TH1F("fCEvents","conter",18,-0.5,17.5);
684  fCEvents->SetStats(kTRUE);
685  fCEvents->GetXaxis()->SetBinLabel(1,"X1");
686  fCEvents->GetXaxis()->SetBinLabel(2,"Analyzed events");
687  fCEvents->GetXaxis()->SetBinLabel(3,"AliAODVertex exists");
688  fCEvents->GetXaxis()->SetBinLabel(4,"TriggerOK");
689  fCEvents->GetXaxis()->SetBinLabel(5,"IsEventSelected");
690  fCEvents->GetXaxis()->SetBinLabel(6,"CascadesHF exists");
691  fCEvents->GetXaxis()->SetBinLabel(7,"MCarray exists");
692  fCEvents->GetXaxis()->SetBinLabel(8,"MCheader exists");
693  fCEvents->GetXaxis()->SetBinLabel(9,"triggerClass!=CINT1");
694  fCEvents->GetXaxis()->SetBinLabel(10,"triggerMask!=kAnyINT");
695  fCEvents->GetXaxis()->SetBinLabel(11,"triggerMask!=kAny");
696  fCEvents->GetXaxis()->SetBinLabel(12,"vtxTitle.Contains(Z)");
697  fCEvents->GetXaxis()->SetBinLabel(13,"vtxTitle.Contains(3D)");
698  fCEvents->GetXaxis()->SetBinLabel(14,"vtxTitle.Doesn'tContain(Z-3D)");
699  fCEvents->GetXaxis()->SetBinLabel(15,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
700  fCEvents->GetXaxis()->SetBinLabel(16,"!IsEventSelected");
701  fCEvents->GetXaxis()->SetBinLabel(17,"triggerMask!=kAnyINT || triggerClass!=CINT1");
702  fCEvents->GetXaxis()->SetBinLabel(18,Form("zVtxMC<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
703  //fCEvents->GetXaxis()->SetTitle("");
704  fCEvents->GetYaxis()->SetTitle("counts");
705 
706  fHTrigger = new TH1F("fHTrigger","counter",18,-0.5,17.5);
707  fHTrigger->SetStats(kTRUE);
708  fHTrigger->GetXaxis()->SetBinLabel(1,"X1");
709  fHTrigger->GetXaxis()->SetBinLabel(2,"kMB");
710  fHTrigger->GetXaxis()->SetBinLabel(3,"kSemiCentral");
711  fHTrigger->GetXaxis()->SetBinLabel(4,"kCentral");
712  fHTrigger->GetXaxis()->SetBinLabel(5,"kINT7");
713  fHTrigger->GetXaxis()->SetBinLabel(6,"kEMC7");
714  //fHTrigger->GetXaxis()->SetBinLabel(7,"Space");
715  fHTrigger->GetXaxis()->SetBinLabel(8,"kMB|kSemiCentral|kCentral");
716  fHTrigger->GetXaxis()->SetBinLabel(9,"kINT7|kEMC7");
717  fHTrigger->GetXaxis()->SetBinLabel(11,"kMB&kSemiCentral");
718  fHTrigger->GetXaxis()->SetBinLabel(12,"kMB&kCentral");
719  fHTrigger->GetXaxis()->SetBinLabel(13,"kINT7&kEMC7");
720 
721  fHCentrality = new TH1F("fHCentrality","conter",100,0.,100.);
722 
723  fOutput->Add(fCEvents);
724  fOutput->Add(fHTrigger);
725  fOutput->Add(fHCentrality);
726 
727  return;
728 }
729 //__________________________________________________________________________
731 {
732  //
733  // Define analyis histograms
734  //
735 
736  //------------------------------------------------
737  // Basic histogram
738  //------------------------------------------------
739  Int_t bins_base[3]= {80 ,20 ,10};
740  Double_t xmin_base[3]={2.286-0.2,0 ,0.00};
741  Double_t xmax_base[3]={2.286+0.2,20. ,100};
742  fHistoLcK0SpMass = new THnSparseF("fHistoLcK0SpMass","",3,bins_base,xmin_base,xmax_base);
744 
745 
746  //------------------------------------------------
747  // checking histograms
748  //------------------------------------------------
749  fHistoBachPt = new TH1F("fHistoBachPt","Bachelor p_{T}",100,0.0,5.0);
750  fOutputAll->Add(fHistoBachPt);
751  fHistod0Bach = new TH1F("fHistod0Bach","Bachelor d_{0}",100,-0.5,0.5);
752  fOutputAll->Add(fHistod0Bach);
753  fHistod0V0 = new TH1F("fHistod0V0","V_{0} d_{0}",100,-0.5,0.5);
754  fOutputAll->Add(fHistod0V0);
755  fHistod0d0 = new TH1F("fHistod0d0","d_{0} d_{0}",100,-0.5,0.5);
756  fOutputAll->Add(fHistod0d0);
757  fHistoV0CosPA=new TH1F("fHistoV0CosPA","V0->Second vertex cospa",100,-1.,1.0);
759  fHistoProbProton=new TH1F("fHistoProbProton","ProbProton",100,0.,1.0);
761  fHistoDecayLength=new TH1F("fHistoDecayLength","Decay Length",100,-0.1,0.1);
763  fHistoK0SMass=new TH1F("fHistoK0SMass","K0S mass",100,0.497-0.05,0.497+0.05);
765 
766  return;
767 }
768 
769 //________________________________________________________________________
770 AliAODRecoCascadeHF* AliAnalysisTaskSELc2pK0sfromAODtracks::MakeCascadeHF(AliAODv0 *v0, AliAODTrack *part, AliAODEvent * aod, AliAODVertex *secVert)
771 {
772  //
773  // Create AliAODRecoCascadeHF object from the argument
774  //
775 
776  if(!v0) return 0x0;
777  if(!part) return 0x0;
778  if(!aod) return 0x0;
779 
780  //------------------------------------------------
781  // PrimaryVertex
782  //------------------------------------------------
783  AliAODVertex *primVertexAOD;
784  Bool_t unsetvtx = kFALSE;
786  primVertexAOD = CallPrimaryVertex(v0,part,aod);
787  if(!primVertexAOD){
788  primVertexAOD = fVtx1;
789  }else{
790  unsetvtx = kTRUE;
791  }
792  }else{
793  primVertexAOD = fVtx1;
794  }
795  if(!primVertexAOD) return 0x0;
796  Double_t posprim[3]; primVertexAOD->GetXYZ(posprim);
797 
798  //------------------------------------------------
799  // DCA between tracks
800  //------------------------------------------------
801  AliESDtrack *esdtrack = new AliESDtrack((AliVTrack*)part);
802 
803  AliNeutralTrackParam *trackV0=NULL;
804  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(v0);
805  if(trackVV0) trackV0 = new AliNeutralTrackParam(trackVV0);
806 
807  Double_t xdummy, ydummy;
808  Double_t dca = esdtrack->GetDCA(trackV0,fBzkG,xdummy,ydummy);
809 
810 
811  //------------------------------------------------
812  // Propagate all tracks to the secondary vertex and calculate momentum there
813  //------------------------------------------------
814 
815  Double_t d0z0bach[2],covd0z0bach[3];
816  if(sqrt(pow(secVert->GetX(),2)+pow(secVert->GetY(),2))<1.){
817  part->PropagateToDCA(secVert,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
818  trackV0->PropagateToDCA(secVert,fBzkG,kVeryBig);
819  }else{
820  part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
821  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig);
822  }
823  Double_t momv0_new[3]={-9999,-9999,-9999.};
824  trackV0->GetPxPyPz(momv0_new);
825 
826  Double_t px[2],py[2],pz[2];
827  px[0] = part->Px(); py[0] = part->Py(); pz[0] = part->Pz();
828  px[1] = momv0_new[0]; py[1] = momv0_new[1]; pz[1] = momv0_new[2];
829 
830  //------------------------------------------------
831  // d0
832  //------------------------------------------------
833  Double_t d0[3],d0err[3];
834 
835  part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0bach,covd0z0bach);
836  d0[0]= d0z0bach[0];
837  d0err[0] = TMath::Sqrt(covd0z0bach[0]);
838 
839  Double_t d0z0v0[2],covd0z0v0[3];
840  trackV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0v0,covd0z0v0);
841  d0[1]= d0z0v0[0];
842  d0err[1] = TMath::Sqrt(covd0z0v0[0]);
843 
844  //------------------------------------------------
845  // Create AliAODRecoCascadeHF
846  //------------------------------------------------
847  Short_t charge = part->Charge();
848  AliAODRecoCascadeHF *theCascade = new AliAODRecoCascadeHF(secVert,charge,px,py,pz,d0,d0err,dca);
849  if(!theCascade)
850  {
851  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
852  if(esdtrack) delete esdtrack;
853  if(trackV0) delete trackV0;
854  return 0x0;
855  }
856  theCascade->SetOwnPrimaryVtx(primVertexAOD);
857  UShort_t id[2]={(UShort_t)part->GetID(),(UShort_t)trackV0->GetID()};
858  theCascade->SetProngIDs(2,id);
859 
860  theCascade->GetSecondaryVtx()->AddDaughter(part);
861  theCascade->GetSecondaryVtx()->AddDaughter(v0);
862 
863  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
864  if(esdtrack) delete esdtrack;
865  if(trackV0) delete trackV0;
866 
867  return theCascade;
868 }
869 
870 //________________________________________________________________________
871 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::CallPrimaryVertex(AliAODv0 *v0, AliAODTrack *trk, AliAODEvent* aod)
872 {
873  //
874  // Make an array of tracks which should not be used in primary vertex calculation and
875  // Call PrimaryVertex function
876  //
877 
878  TObjArray *TrackArray = new TObjArray(3);
879 
880  AliESDtrack *cptrk1 = new AliESDtrack((AliVTrack*)trk);
881  TrackArray->AddAt(cptrk1,0);
882 
883  AliESDtrack *cascptrack = new AliESDtrack((AliVTrack*)v0->GetDaughter(0));
884  TrackArray->AddAt(cascptrack,1);
885  AliESDtrack *cascntrack = new AliESDtrack((AliVTrack*)v0->GetDaughter(1));
886  TrackArray->AddAt(cascntrack,2);
887 
888  AliAODVertex *newvert = PrimaryVertex(TrackArray,aod);
889 
890  for(Int_t i=0;i<3;i++)
891  {
892  AliESDtrack *tesd = (AliESDtrack*)TrackArray->UncheckedAt(i);
893  delete tesd;
894  }
895  TrackArray->Clear();
896  delete TrackArray;
897 
898  return newvert;
899 }
900 
901 //________________________________________________________________________
902 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::PrimaryVertex(const TObjArray *trkArray,
903  AliVEvent *event)
904 {
905  //
906  //Used only for pp
907  //copied from AliAnalysisVertexingHF (except for the following 3 lines)
908  //
909 
910  Bool_t fRecoPrimVtxSkippingTrks = kTRUE;
911  Bool_t fRmTrksFromPrimVtx = kFALSE;
912 
913  AliESDVertex *vertexESD = 0;
914  AliAODVertex *vertexAOD = 0;
915 
916  //vertexESD = new AliESDVertex(*fV1);
917 
918 
919  if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
920  // primary vertex from the input event
921 
922  vertexESD = new AliESDVertex(*fV1);
923 
924  } else {
925  // primary vertex specific to this candidate
926 
927  Int_t nTrks = trkArray->GetEntriesFast();
928  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
929 
930  if(fRecoPrimVtxSkippingTrks) {
931  // recalculating the vertex
932 
933  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
934  Float_t diamondcovxy[3];
935  event->GetDiamondCovXY(diamondcovxy);
936  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
937  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
938  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
939  vertexer->SetVtxStart(diamond);
940  delete diamond; diamond=NULL;
941  if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
942  vertexer->SetOnlyFitter();
943  }
944  Int_t skipped[1000];
945  Int_t nTrksToSkip=0,id;
946  AliExternalTrackParam *t = 0;
947  for(Int_t i=0; i<nTrks; i++) {
948  t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
949  id = (Int_t)t->GetID();
950  if(id<0) continue;
951  skipped[nTrksToSkip++] = id;
952  }
953  // TEMPORARY FIX
954  // For AOD, skip also tracks without covariance matrix
955  Double_t covtest[21];
956  for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
957  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
958  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
959  id = (Int_t)vtrack->GetID();
960  if(id<0) continue;
961  skipped[nTrksToSkip++] = id;
962  }
963  }
964  for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
965  //
966  vertexer->SetSkipTracks(nTrksToSkip,skipped);
967  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
968 
969  } else if(fRmTrksFromPrimVtx && nTrks>0) {
970  // removing the prongs tracks
971 
972  TObjArray rmArray(nTrks);
973  UShort_t *rmId = new UShort_t[nTrks];
974  AliESDtrack *esdTrack = 0;
975  AliESDtrack *t = 0;
976  for(Int_t i=0; i<nTrks; i++) {
977  t = (AliESDtrack*)trkArray->UncheckedAt(i);
978  esdTrack = new AliESDtrack(*t);
979  rmArray.AddLast(esdTrack);
980  if(esdTrack->GetID()>=0) {
981  rmId[i]=(UShort_t)esdTrack->GetID();
982  } else {
983  rmId[i]=9999;
984  }
985  }
986  Float_t diamondxy[2]={static_cast<Float_t>(event->GetDiamondX()),static_cast<Float_t>(event->GetDiamondY())};
987  vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
988  delete [] rmId; rmId=NULL;
989  rmArray.Delete();
990 
991  }
992 
993  delete vertexer; vertexer=NULL;
994  if(!vertexESD) return vertexAOD;
995  if(vertexESD->GetNContributors()<=0) {
996  //AliDebug(2,"vertexing failed");
997  delete vertexESD; vertexESD=NULL;
998  return vertexAOD;
999  }
1000 
1001 
1002  }
1003 
1004  // convert to AliAODVertex
1005  Double_t pos[3],cov[6],chi2perNDF;
1006  vertexESD->GetXYZ(pos); // position
1007  vertexESD->GetCovMatrix(cov); //covariance matrix
1008  chi2perNDF = vertexESD->GetChi2toNDF();
1009  delete vertexESD; vertexESD=NULL;
1010 
1011  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1012 
1013  return vertexAOD;
1014 }
1015 
1016 //________________________________________________________________________
1017 AliAODVertex* AliAnalysisTaskSELc2pK0sfromAODtracks::ReconstructSecondaryVertex(AliAODv0 *v0, AliAODTrack *part, AliAODEvent * aod)
1018 {
1019  //
1020  // Reconstruct secondary vertex from trkArray (Copied from AliAnalysisVertexingHF)
1021  //
1022  //------------------------------------------------
1023  // PrimaryVertex
1024  //------------------------------------------------
1025  AliAODVertex *primVertexAOD;
1026  Bool_t unsetvtx = kFALSE;
1028  primVertexAOD = CallPrimaryVertex(v0,part,aod);
1029  if(!primVertexAOD){
1030  primVertexAOD = fVtx1;
1031  }else{
1032  unsetvtx = kTRUE;
1033  }
1034  }else{
1035  primVertexAOD = fVtx1;
1036  }
1037  if(!primVertexAOD) return 0x0;
1038 
1039  //------------------------------------------------
1040  // Secondary vertex
1041  //------------------------------------------------
1042 
1043  Double_t LcPx = part->Px()+v0->Px();
1044  Double_t LcPy = part->Py()+v0->Py();
1045  Double_t LcPt = TMath::Sqrt(LcPx*LcPx+LcPy*LcPy);
1046 
1047  Double_t d0z0[2],covd0z0[3];
1048  part->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
1049  Double_t x0 = primVertexAOD->GetX();
1050  Double_t y0 = primVertexAOD->GetY();
1051  Double_t px0 = LcPx/LcPt;
1052  Double_t py0 = LcPy/LcPt;
1053  Double_t tx[3];
1054  part->GetXYZ(tx);
1055  Double_t x1 = tx[0];
1056  Double_t y1 = tx[1];
1057  part->GetPxPyPz(tx);
1058  Double_t px1 = tx[0];
1059  Double_t py1 = tx[1];
1060  Double_t pt1 = sqrt(px1*px1+py1*py1);
1061  px1 = px1/pt1;
1062  py1 = py1/pt1;
1063 
1064  Double_t dx = x0 - x1;
1065  Double_t dy = y0 - y1;
1066 
1067  Double_t Delta = -px0*py1+py0*px1;
1068  Double_t a0 = -9999.;
1069  if(Delta!=0)
1070  {
1071  a0 = 1./Delta * (py1 * dx - px1 * dy);
1072  }
1073  Double_t neovertx = x0 + a0 * px0;
1074  Double_t neoverty = y0 + a0 * py0;
1075  Double_t z0 = primVertexAOD->GetZ();
1076  Double_t neovertz = z0 + TMath::Abs(a0)*part->Pz()/part->Pt();
1077 
1078  if(unsetvtx) delete primVertexAOD; primVertexAOD=NULL;
1079 
1080  Double_t pos[3],cov[6],chi2perNDF;
1081  pos[0]=neovertx;
1082  pos[1]=neoverty;
1083  pos[2]=neovertz;
1084  cov[0]=0.0;
1085  cov[1]=0.0;
1086  cov[2]=0.0;
1087  cov[3]=0.0;
1088  cov[4]=0.0;
1089  cov[5]=0.0;
1090  chi2perNDF=0.0;
1091  AliAODVertex *secVert = new AliAODVertex(pos,cov,chi2perNDF);
1092  if(!secVert){
1093  return 0x0;
1094  }
1095  return secVert;
1096 }
Int_t charge
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Bool_t GetUseCombined()
Definition: AliAODPidHF.h:165
AliAODVertex * PrimaryVertex(const TObjArray *trkArray, AliVEvent *event)
TH1F * fHistoProbProton
! Probability to be proton histogram
AliAODv0 * Getv0() const
TTree * fVariablesTree
flag to decide whether to write the candidate variables on a tree variables
void FillROOTObjects(AliAODRecoCascadeHF *lcobj, AliAODMCParticle *mcpart, AliAODMCParticle *mcdau1, AliAODMCParticle *mcdau2, Int_t mcnused)
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:241
Double_t InvMassLctoK0sP() const
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:231
TH1F * fHCentrality
! Histogram to check Centrality
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
void SetProngIDs(Int_t nIDs, UShort_t *id)
THnSparse * fHistoLcK0SpMass
Stores trigger information.
AliPIDResponse * GetPidResponse() const
Definition: AliAODPidHF.h:160
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)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Double_t DecayLengthXY() const
Bool_t IsEventSelected(AliVEvent *event)
Bool_t fIsMB
Reconstruct primary vertex excluding candidate tracks.
TList * fOutputAll
! User Output slot 3 //analysis histograms
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:253
TH1F * fCEvents
! Histogram to check selected events
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:295
AliAODRecoCascadeHF * MakeCascadeHF(AliAODv0 *casc, AliAODTrack *trk, AliAODEvent *aod, AliAODVertex *vert)
AliAODVertex * CallPrimaryVertex(AliAODv0 *v0, AliAODTrack *trk, AliAODEvent *evt)
Bool_t fIsEventSelected
Cuts - sent to output slot 2.
AliAODVertex * ReconstructSecondaryVertex(AliAODv0 *casc, AliAODTrack *trk, AliAODEvent *aod)
TH1F * fHistoV0CosPA
! V0 cosine pointing angle to primary vertex