Fix for leak
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskCombinHF.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2018, 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  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: $ */
17
18 //*************************************************************************
19 // Class AliAnalysisTaskCombinHF
20 // AliAnalysisTaskSE to build D meson candidates by combining tracks
21 //  background is computed LS and track rotations is
22 // Authors: F. Prino, A. Rossi
23 /////////////////////////////////////////////////////////////
24
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TH3F.h>
29 #include <THnSparse.h>
30
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
33 #include "AliPIDResponse.h"
34 #include "AliAODHandler.h"
35 #include "AliAODEvent.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODVertex.h"
39 #include "AliAODTrack.h"
40 #include "AliVertexingHFUtils.h"
41 #include "AliAnalysisTaskCombinHF.h"
42
43 ClassImp(AliAnalysisTaskCombinHF)
44
45
46 //________________________________________________________________________
47 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
48   AliAnalysisTaskSE(),
49   fOutput(0x0), 
50   fHistNEvents(0x0),
51   fHistTrackStatus(0x0),
52   fHistCheckOrigin(0x0),
53   fHistCheckOriginSel(0x0),
54   fHistCheckDecChan(0x0),
55   fHistCheckDecChanAcc(0x0),
56   fPtVsYGen(0x0),
57   fPtVsYGenLimAcc(0x0),
58   fPtVsYGenAcc(0x0),
59   fPtVsYReco(0x0),
60   fMassVsPtVsY(0x0),
61   fMassVsPtVsYRot(0x0),
62   fMassVsPtVsYLSpp(0x0),
63   fMassVsPtVsYLSmm(0x0),
64   fMassVsPtVsYSig(0x0),
65   fMassVsPtVsYRefl(0x0),
66   fMassVsPtVsYBkg(0x0),
67   fNSelected(0x0),
68   fNormRotated(0x0),
69   fDeltaMass(0x0),
70   fDeltaMassFullAnalysis(0x0),
71   fFilterMask(BIT(4)),
72   fTrackCutsAll(0x0),
73   fTrackCutsPion(0x0),
74   fTrackCutsKaon(0x0),
75   fPidHF(new AliAODPidHF()),
76   fAnalysisCuts(0x0),
77   fMinMass(1.750),
78   fMaxMass(2.150),
79   fEtaAccCut(0.9),
80   fPtAccCut(0.1),
81   fNRotations(9),
82   fMinAngleForRot(5*TMath::Pi()/6),
83   fMaxAngleForRot(7*TMath::Pi()/6),
84   fMinAngleForRot3(2*TMath::Pi()/6),
85   fMaxAngleForRot3(4*TMath::Pi()/6),
86   fCounter(0x0),
87   fMeson(kDzero),
88   fReadMC(kFALSE),
89   fPromptFeeddown(kPrompt),
90   fGoUpToQuark(kTRUE),
91   fFullAnalysis(0)  
92 {
93   // default constructor
94 }
95
96 //________________________________________________________________________
97 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts):
98   AliAnalysisTaskSE("DmesonCombin"),
99   fOutput(0x0), 
100   fHistNEvents(0x0),
101   fHistTrackStatus(0x0),
102   fHistCheckOrigin(0x0),
103   fHistCheckOriginSel(0x0),
104   fHistCheckDecChan(0x0),
105   fHistCheckDecChanAcc(0x0),
106   fPtVsYGen(0x0),
107   fPtVsYGenLimAcc(0x0),
108   fPtVsYGenAcc(0x0),
109   fPtVsYReco(0x0),
110   fMassVsPtVsY(0x0),
111   fMassVsPtVsYRot(0x0),
112   fMassVsPtVsYLSpp(0x0),
113   fMassVsPtVsYLSmm(0x0),
114   fMassVsPtVsYSig(0x0),
115   fMassVsPtVsYRefl(0x0),
116   fMassVsPtVsYBkg(0x0),
117   fNSelected(0x0),
118   fNormRotated(0x0),
119   fDeltaMass(0x0),
120   fDeltaMassFullAnalysis(0x0),
121   fFilterMask(BIT(4)),
122   fTrackCutsAll(0x0),
123   fTrackCutsPion(0x0),
124   fTrackCutsKaon(0x0),
125   fPidHF(new AliAODPidHF()),
126   fAnalysisCuts(analysiscuts),
127   fMinMass(1.750),
128   fMaxMass(2.150),
129   fEtaAccCut(0.9),
130   fPtAccCut(0.1),
131   fNRotations(9),
132   fMinAngleForRot(5*TMath::Pi()/6),
133   fMaxAngleForRot(7*TMath::Pi()/6),
134   fMinAngleForRot3(2*TMath::Pi()/6),
135   fMaxAngleForRot3(4*TMath::Pi()/6),
136   fCounter(0x0),
137   fMeson(meson),
138   fReadMC(kFALSE),
139   fPromptFeeddown(1),
140   fGoUpToQuark(kTRUE),
141   fFullAnalysis(0)
142
143 {
144   // standard constructor
145
146   DefineOutput(1,TList::Class());  //My private output
147   DefineOutput(2,AliNormalizationCounter::Class());
148  }
149
150 //________________________________________________________________________
151 AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
152 {
153   //
154   // Destructor
155   //
156   delete fOutput;
157   delete fHistNEvents;
158   delete fHistTrackStatus;
159   delete fHistCheckOrigin;
160   delete fHistCheckOriginSel;
161   delete fHistCheckDecChan;
162   delete fHistCheckDecChanAcc;
163   delete fPtVsYGen;
164   delete fPtVsYGenLimAcc;
165   delete fPtVsYGenAcc;
166   delete fPtVsYReco;
167   delete fMassVsPtVsY; 
168   delete fMassVsPtVsYLSpp;
169   delete fMassVsPtVsYLSmm;
170   delete fMassVsPtVsYRot;
171   delete fMassVsPtVsYSig;
172   delete fMassVsPtVsYRefl;
173   delete fMassVsPtVsYBkg;
174   delete fNSelected;
175   delete fNormRotated;
176   delete fDeltaMass;
177   delete fCounter;
178   delete fTrackCutsAll;
179   delete fTrackCutsPion;
180   delete fTrackCutsKaon;
181   delete fPidHF;
182   delete fAnalysisCuts;
183 }
184
185 //________________________________________________________________________
186 void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
187 {
188   // Create the output container
189   //
190   if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
191
192   fOutput = new TList();
193   fOutput->SetOwner();
194   fOutput->SetName("OutputHistos");
195
196   fHistNEvents = new TH1F("hNEvents", "number of events ",8,-0.5,7.5);
197   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
198   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
199   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
200   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
201   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
202   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
203   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
204   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
205
206   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
207   fHistNEvents->Sumw2();
208   fHistNEvents->SetMinimum(0);
209   fOutput->Add(fHistNEvents);
210
211   fHistTrackStatus  = new TH1F("hTrackStatus", "",8,-0.5,7.5);
212   fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
213   fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
214   fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
215   fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
216   fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
217   fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
218   fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
219   fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
220   fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
221   fHistTrackStatus->Sumw2();
222   fHistTrackStatus->SetMinimum(0);
223   fOutput->Add(fHistTrackStatus);
224
225   if(fReadMC){
226
227     fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
228     fHistCheckOrigin->Sumw2();
229     fHistCheckOrigin->SetMinimum(0);
230     fOutput->Add(fHistCheckOrigin);
231
232     fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
233     fHistCheckOriginSel->Sumw2();
234     fHistCheckOriginSel->SetMinimum(0);
235     fOutput->Add(fHistCheckOriginSel);
236
237     fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
238     fHistCheckDecChan->Sumw2();
239     fHistCheckDecChan->SetMinimum(0);
240     fOutput->Add(fHistCheckDecChan);
241
242     fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
243     fHistCheckDecChanAcc->Sumw2();
244     fHistCheckDecChanAcc->SetMinimum(0);
245     fOutput->Add(fHistCheckDecChanAcc);
246
247     fPtVsYGen= new TH2F("hPtVsYGen","",20,0.,10.,20,-1.,1.);
248     fPtVsYGen->Sumw2();
249     fPtVsYGen->SetMinimum(0);
250     fOutput->Add(fPtVsYGen);
251
252     fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",20,0.,10.,20,-1.,1.);
253     fPtVsYGenLimAcc->Sumw2();
254     fPtVsYGenLimAcc->SetMinimum(0);
255     fOutput->Add(fPtVsYGenLimAcc);
256
257     fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",20,0.,10.,20,-1.,1.);
258     fPtVsYGenAcc->Sumw2();
259     fPtVsYGenAcc->SetMinimum(0);
260     fOutput->Add(fPtVsYGenAcc);
261
262     fPtVsYReco= new TH2F("hPtVsYReco","",20,0.,10.,20,-1.,1.);
263     fPtVsYReco->Sumw2();
264     fPtVsYReco->SetMinimum(0);
265     fOutput->Add(fPtVsYReco);
266   }
267
268
269  Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.;
270   Double_t maxm=fMinMass+nMassBins*0.001;
271   fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
272   fMassVsPtVsY->Sumw2();
273   fMassVsPtVsY->SetMinimum(0);
274   fOutput->Add(fMassVsPtVsY);
275
276   fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
277   fMassVsPtVsYRot->Sumw2();
278   fMassVsPtVsYRot->SetMinimum(0);
279   fOutput->Add(fMassVsPtVsYRot);
280
281   if(fMeson==kDzero){
282     fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
283     fMassVsPtVsYLSpp->Sumw2();
284     fMassVsPtVsYLSpp->SetMinimum(0);
285     fOutput->Add(fMassVsPtVsYLSpp);
286     fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
287     fMassVsPtVsYLSmm->Sumw2();
288     fMassVsPtVsYLSmm->SetMinimum(0);
289     fOutput->Add(fMassVsPtVsYLSmm);
290   }
291
292   fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
293   fMassVsPtVsYSig->Sumw2();
294   fMassVsPtVsYSig->SetMinimum(0);
295   fOutput->Add(fMassVsPtVsYSig);
296
297   fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
298   fMassVsPtVsYRefl->Sumw2();
299   fMassVsPtVsYRefl->SetMinimum(0);
300   fOutput->Add(fMassVsPtVsYRefl);
301
302   fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
303   fMassVsPtVsYBkg->Sumw2();
304   fMassVsPtVsYBkg->SetMinimum(0);
305   fOutput->Add(fMassVsPtVsYBkg);
306
307   fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
308   fNSelected->Sumw2();
309   fNSelected->SetMinimum(0);
310   fOutput->Add(fNSelected);
311
312   fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
313   fNormRotated->Sumw2();
314   fNormRotated->SetMinimum(0);
315   fOutput->Add(fNormRotated);
316
317   fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
318   fDeltaMass->Sumw2();
319   fDeltaMass->SetMinimum(0);
320   fOutput->Add(fDeltaMass);
321
322   Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
323   Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
324   Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
325   fDeltaMassFullAnalysis=new THnSparseF("fDeltaMassFullAnalysis","fDeltaMassFullAnalysis;inv mass (GeV/c);#Delta inv mass (GeV/c) ; p_{T}^{D} (GeV/c); #Delta p_{T} (GeV/c); daughter angle (2prongs) (rad);",5,binSparseDMassRot,edgeLowSparseDMassRot,edgeHighSparseDMassRot);
326   fOutput->Add(fDeltaMassFullAnalysis);
327
328   //Counter for Normalization
329   fCounter = new AliNormalizationCounter("NormalizationCounter");
330   fCounter->Init();
331
332   PostData(1,fOutput); 
333   PostData(2,fCounter);   
334 }
335
336 //________________________________________________________________________
337 void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
338   //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
339
340   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
341   if(!aod && AODEvent() && IsStandardAOD()) {
342     // In case there is an AOD handler writing a standard AOD, use the AOD 
343     // event in memory rather than the input (ESD) event.    
344     aod = dynamic_cast<AliAODEvent*> (AODEvent());
345   }
346   if(!aod){
347     printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
348     return;
349   }
350
351   // fix for temporary bug in ESDfilter 
352   // the AODs with null vertex pointer didn't pass the PhysSel
353   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
354   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
355   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
356   AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
357   fPidHF->SetPidResponse(pidResp);
358   
359  
360   fHistNEvents->Fill(0); // count event
361   // Post the data already here
362   PostData(1,fOutput);
363
364   fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
365  
366   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
367   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
368   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
369   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
370   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
371   //  if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6);
372   if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7); 
373
374   if(!isEvSel)return;
375  
376   fHistNEvents->Fill(1);
377   
378   TClonesArray *arrayMC=0;
379   AliAODMCHeader *mcHeader=0;
380   if(fReadMC){
381     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
382     if(!arrayMC) {
383       printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
384       return;
385     }
386     
387     // load MC header
388     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
389     if(!mcHeader) {
390       printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
391       return;
392     }
393     FillGenHistos(arrayMC);
394   }
395
396
397   Int_t ntracks=aod->GetNTracks();
398
399   // select and flag tracks
400   UChar_t* status = new UChar_t[ntracks];
401   for(Int_t iTr=0; iTr<ntracks; iTr++){
402     status[iTr]=0;
403     AliAODTrack* track=aod->GetTrack(iTr);
404     if(IsTrackSelected(track)) status[iTr]+=1;
405     if(IsKaon(track)) status[iTr]+=2;
406     if(IsPion(track)) status[iTr]+=4;
407     fHistTrackStatus->Fill(status[iTr]);
408   }
409
410   // build the combinatorics
411   Int_t nSelected=0;
412   Int_t nFiltered=0;
413   Double_t dummypos[3]={0.,0.,0.};
414   AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
415   AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
416   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
417   Double_t d02[2]={0.,0.};
418   Double_t d03[3]={0.,0.,0.};
419   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
420   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
421   UInt_t pdg0[2]={321,211};
422   UInt_t pdgp[3]={321,211,211};
423   //  UInt_t pdgs[3]={321,321,211};
424   Double_t tmpp[3];
425   Double_t px[3],py[3],pz[3];
426   Int_t dgLabels[3];
427
428   for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
429     AliAODTrack* trK=aod->GetTrack(iTr1);
430     if((status[iTr1] & 1)==0) continue;
431     if((status[iTr1] & 2)==0) continue;
432     Int_t chargeK=trK->Charge();
433     trK->GetPxPyPz(tmpp);
434     px[0] = tmpp[0]; 
435     py[0] = tmpp[1]; 
436     pz[0] = tmpp[2]; 
437     dgLabels[0]=trK->GetLabel();
438     for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
439       if((status[iTr2] & 1)==0) continue;
440       if((status[iTr2] & 4)==0) continue;
441       if(iTr1==iTr2) continue;
442       AliAODTrack* trPi1=aod->GetTrack(iTr2);
443       Int_t chargePi1=trPi1->Charge();
444       trPi1->GetPxPyPz(tmpp);
445       px[1] = tmpp[0]; 
446       py[1] = tmpp[1]; 
447       pz[1] = tmpp[2]; 
448       dgLabels[1]=trPi1->GetLabel();
449       if(chargePi1==chargeK){
450         if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
451         continue;
452       }
453       if(fMeson==kDzero){
454         nFiltered++;
455         v2->AddDaughter(trK);
456         v2->AddDaughter(trPi1);
457         tmpRD2->SetSecondaryVtx(v2);
458         Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
459         v2->RemoveDaughters();
460         if(ok) nSelected++;
461       }else{
462         for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
463           if((status[iTr3] & 1)==0) continue;
464           if((status[iTr3] & 4)==0) continue;
465           if(iTr1==iTr3) continue;
466           AliAODTrack* trPi2=aod->GetTrack(iTr3);
467           Int_t chargePi2=trPi2->Charge();
468           if(chargePi2==chargeK) continue;
469           nFiltered++;
470           trPi2->GetPxPyPz(tmpp);
471           px[2] = tmpp[0]; 
472           py[2] = tmpp[1]; 
473           pz[2] = tmpp[2]; 
474           dgLabels[2]=trPi2->GetLabel();
475           v3->AddDaughter(trK);
476           v3->AddDaughter(trPi1);
477           v3->AddDaughter(trPi2);
478           tmpRD3->SetSecondaryVtx(v3);
479           Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
480           v3->RemoveDaughters();
481           if(ok) nSelected++;
482         }
483       }
484     }
485   }
486
487   delete [] status;
488   delete v2;
489   delete v3;
490   delete tmpRD2;
491   delete tmpRD3;
492
493   fNSelected->Fill(nSelected);
494
495   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
496   fCounter->StoreCandidates(aod,nSelected,kFALSE);
497
498   PostData(1,fOutput); 
499   PostData(2,fCounter);    
500
501   return;
502 }
503
504 //________________________________________________________________________
505 void AliAnalysisTaskCombinHF::FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){
506   // Fill histos for LS candidates
507     
508   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
509   Double_t pt = tmpRD->Pt();
510   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
511   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
512     Double_t rapid = tmpRD->Y(pdgD);
513     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
514       if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
515       else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
516     }
517   }
518   return;
519 }
520
521 //________________________________________________________________________
522 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
523   // Fill histos with generated quantities
524   Int_t totPart=arrayMC->GetEntriesFast();
525   Int_t thePDG=411;
526   Int_t nProng=3;
527   if(fMeson==kDzero){
528     thePDG=421;
529     nProng=2;
530   }
531   for(Int_t ip=0; ip<totPart; ip++){
532     AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
533     if(TMath::Abs(part->GetPdgCode())==thePDG){
534       Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
535       fHistCheckOrigin->Fill(orig);
536       if(fPromptFeeddown==kFeeddown && orig!=5) continue;
537       else if(fPromptFeeddown==kPrompt && orig!=4) continue;
538       else if(fPromptFeeddown==kBoth && orig<4) continue;
539       fHistCheckOriginSel->Fill(orig);
540       Int_t deca=0;
541       Bool_t isGoodDecay=kFALSE;
542       Int_t labDau[4]={-1,-1,-1,-1};
543       if(fMeson==kDzero){
544         deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
545         if(part->GetNDaughters()!=2) continue;
546         if(deca==1) isGoodDecay=kTRUE;
547       }else if(fMeson==kDplus){ 
548         deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
549         if(deca>0) isGoodDecay=kTRUE;
550       }
551       fHistCheckDecChan->Fill(deca);
552       if(labDau[0]==-1){
553         //      printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
554         continue; //protection against unfilled array of labels
555       }
556       Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
557       if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
558       if(isGoodDecay){
559         Double_t ptgen=part->Pt();
560         Double_t ygen=part->Y();
561         if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
562           fPtVsYGen->Fill(ptgen,ygen);
563           if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen);
564           if(isInAcc) fPtVsYGenAcc->Fill(ptgen,ygen);
565         }
566       }
567     }
568   }
569 }
570
571 //________________________________________________________________________
572 Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels){
573   // Fill histos for candidates with proper charge sign
574
575   Bool_t accept=kFALSE;
576
577   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
578   Double_t pt = tmpRD->Pt();
579   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
580   Double_t mass=TMath::Sqrt(minv2);
581
582   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
583     Double_t rapid = tmpRD->Y(pdgD);
584     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
585       fMassVsPtVsY->Fill(mass,pt,rapid);
586       accept=kTRUE;
587       if(fReadMC){
588         Int_t signPdg[3]={0,0,0};
589         for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
590         Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
591         if(labD>=0){
592           AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
593           if(part){
594             Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
595             if(pdgCode==321){
596               AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
597               if(dmes){
598                 Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,dmes,fGoUpToQuark);
599                 if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) {
600                   fPtVsYReco->Fill(dmes->Pt(),dmes->Y());
601                 }
602               }
603               fMassVsPtVsYSig->Fill(mass,pt,rapid);
604             }else{
605               fMassVsPtVsYRefl->Fill(mass,pt,rapid);
606             }
607           }
608         }else{
609           fMassVsPtVsYBkg->Fill(mass,pt,rapid);
610         }
611       }
612     }
613   }
614
615   Int_t nRotated=0;
616   Double_t massRot=0;// calculated later only if candidate is acceptable
617   Double_t angleProngXY;
618   if(TMath::Abs(pdgD)==421)angleProngXY=TMath::ACos((px[0]*px[1]+py[0]*py[1])/TMath::Sqrt((px[0]*px[0]+py[0]*py[0])*(px[1]*px[1]+py[1]*py[1])));
619   else {
620     angleProngXY=TMath::ACos(((px[0]+px[1])*px[2]+(py[0]+py[1])*py[2])/TMath::Sqrt(((px[0]+px[1])*(px[0]+px[1])+(py[0]+py[1])*(py[0]+py[1]))*(px[2]*px[2]+py[2]*py[2])));
621   }
622   Double_t ptOrig=pt;
623
624
625   Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
626   Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
627
628   for(Int_t irot=0; irot<fNRotations; irot++){
629     Double_t phirot=fMinAngleForRot+rotStep*irot;
630     Double_t tmpx=px[0];
631     Double_t tmpy=py[0];
632     Double_t tmpx2=px[2];
633     Double_t tmpy2=py[2];
634     px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
635     py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
636     if(pdgD==411 || pdgD==431){
637       Double_t phirot2=fMinAngleForRot3+rotStep3*irot;
638       px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
639       py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
640     }
641     tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
642     pt = tmpRD->Pt();
643     minv2 = tmpRD->InvMass2(nProngs,pdgdau);
644     if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
645       Double_t rapid = tmpRD->Y(pdgD);
646       if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
647         massRot=TMath::Sqrt(minv2);
648         fMassVsPtVsYRot->Fill(massRot,pt,rapid);
649         nRotated++;
650         fDeltaMass->Fill(massRot-mass);
651         if(fFullAnalysis){
652           Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
653           fDeltaMassFullAnalysis->Fill(pointRot);
654         }
655       }
656     }
657     px[0]=tmpx;
658     py[0]=tmpy;
659     if(pdgD==411 || pdgD==431){
660       px[2]=tmpx2;
661       py[2]=tmpy2;
662     }
663   }
664   fNormRotated->Fill(nRotated);
665
666   return accept;
667
668 }
669 //________________________________________________________________________
670 Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
671   // track selection cuts
672
673   if(track->Charge()==0) return kFALSE;
674   if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
675   if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;    
676   return kTRUE;
677 }
678 //________________________________________________________________________
679 Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){
680   // kaon selection cuts
681
682   if(!fPidHF) return kTRUE;
683   Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);  
684   if(isKaon>=0 && SelectAODTrack(track,fTrackCutsKaon)) return kTRUE;
685   return kFALSE;
686 }
687 //________________________________________________________________________
688 Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){
689   // pion selection cuts
690
691   if(!fPidHF) return kTRUE;
692   Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
693   if(isPion>=0&& SelectAODTrack(track,fTrackCutsPion)) return kTRUE;
694   return kFALSE;
695 }
696 //________________________________________________________________________
697 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
698   // AOD track selection
699
700   if(!cuts) return kTRUE;
701
702   AliESDtrack esdTrack(track);
703   // set the TPC cluster info
704   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
705   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
706   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
707   if(!cuts->IsSelected(&esdTrack)) return kFALSE; 
708
709   return kTRUE;  
710 }
711
712 //_________________________________________________________________
713 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
714   // check if the decay products are in the good eta and pt range
715   for (Int_t iProng = 0; iProng<nProng; iProng++){
716     AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
717     if(!mcPartDaughter) return kFALSE;
718     Double_t eta = mcPartDaughter->Eta();
719     Double_t pt = mcPartDaughter->Pt();
720     if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
721   }
722   return kTRUE;
723 }
724
725 //_________________________________________________________________
726 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
727 {
728   // Terminate analysis
729   //
730   if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
731   fOutput = dynamic_cast<TList*> (GetOutputData(1));
732   if (!fOutput) {     
733     printf("ERROR: fOutput not available\n");
734     return;
735   }
736   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
737   if(fHistNEvents){
738     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
739   }else{
740     printf("ERROR: fHistNEvents not available\n");
741     return;
742   }
743   return;
744 }
745