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